Open Case Studies : Exploring global patterns of dietary behaviors associated with health risk

Disclaimer: The purpose of the Open Case Studies project is to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given dataset, and should not be used in the context of making policy decisions without external consultation from scientific experts.

Motivation

An article was recently published in the lancet journal that evaluates global dietary trends and the relationship of these dietary factors with mortality and fertility.

GBD 2017 Diet Collaborators. Health effects of dietary risks in 195 countries, 1990–2017: a systematic analysis for the Global Burden of Disease Study 2017. The Lancet 393, 1958–1972 (2019).

This article evaluated food consumption patterns in 195 countries for 15 different dietary risk factors that have probable associations with non-communicable disease (NCD). For example, overconsumption of sodium is associated with high blood pressure. These consumption levels were then used to estimate levels of mortality and morbidity due NCD, as well as disability-adjusted life-years (DALYs) attributible to suboptimal consumption of foods related to these dietary risk factors. The authors found that:

“High intake of sodium …, low intake of whole grains …, and low intake of fruits … were the leading dietary risk factors for deaths and DALYs globally and in many countries.”

This figure from the supplement shows the ranking of the 15 dietary risk factors based on the estimated number of attributable deaths and illustartes how the top 3 risk factors are often issues for many different countries.

This case study will evaluate the data reported in this article to explore regional, age, and gender specific differences in dietary consumption patterns around the world in 2017.

Main Questions

Our main questions are:

1. What are the global trends for potentially harmful diets?
2. How do males and females compare?
3. How do different age groups compare for these dietary factors?
4. How do different countries compare? In particular, how does the US compare to other contries in terms of diet trends?

Learning Objectives

In this case study, we’ll walk you through importing data from a pdf, cleaning data, wrangling data, visualizing the data, and comparing two or more groups using well-established and commonly used packages, including stringr, tidyr, dplyr, purrr, and ggplot2. We will especially focus on using packages and functions from the Tidyverse. The Tidyverse is a library of packages created by the chief scientist at RStudio, Hadley Wickham. While some students may be familiar with previous R programming packages, these packages make data science in R especially efficient.

We will begin by loading the packages that we will need:

library(here)
library(readr)
library(dplyr)
library(skimr)
library(pdftools)
library(stringr)
library(magrittr)
library(purrr)
library(tibble)
library(tidyr)
library(ggplot2)
library(ggpubr)
library(forcats)
library(ggrepel)
library(cowplot)

Package	Use
here	to easily load and save data
readr	to import the csv file data
dplyr	to arrange/filter/select/compare specific subsets of the data
skimr	to get an overview of data
pdftools	to read a pdf into R
stringr	to manipulate the text within the pdf of the data
magrittr	to use the %<>% pipping operator
purrr	to perform functions on all columns of a tibble
tibble	to create data objects that we can manipulate with dplyr/stringr/tidyr/purrr
tidyr	to separate data within a column into multiple columns
ggplot2	to make visualizations with multiple layers
ggpubr	to easily add regression line equations to plots
forcats	to change details about factors (categorical variables)
ggrepel	to allow labels in figures not to overlap
cowplot	to allow plots to be combined

---

glue | to paste or combine character strings and data together

The first time we use a function, we will use the :: to indicate which package we are using. Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.

Context

Here is an excerpt from the article itself about the context of the work:

Many dietary factors have well-estabished associations with health risk. The authors that generated this dataset identified 15 dietary factors that have probably health risk based on literature search.

Here you can see a table of the sources for the health risks associated with the dietary factors. RCT stands for randomized control trials.

In the article the authors found that most of the mortality associated with each factor is related to cardiovascular disease.

Limitations

There are some important limitations regarding the data from this article to keep in mind. The definition of certain dietary factors varied across some of the collection sources. Intakes of certain healthy foods like vegitables and fruits are likely positively correlated and likely negatively correlated with intakes of unhealthy foods. Much of the data was collected with 24 hour recall surveys which are prone to issues due to inaccuracy of memory recall or other biases such as a tendency for some people to report healthier behaviriors. The guidelines in the table are based are not parsed by gender even though it is known that there are different dietary requirements for optimal health for certain nutrients. The article discusses some limitations about accounting for overall food consumption when calculting consumption of particular foods:

"To remove the effect of energy intake as a potential confounder and address measurement error in dietary assessment tools, most cohorts have adjusted for total energy intake in their statistical models. This energy adjustment means that diet components are defined as risks in terms of the share of diet and not as absolute levels of exposure. In other words, an increase in intake of foods and macronutrients should be compensated by a decrease in intake of other dietary factors to hold total energy intake constant. Thus, the relative risk of change in each component of diet depends on the other components for which it is substituted. However, the relative risks estimated from meta-analyses of cohort studies do not generally specify the type of substitution.

There are also important nuances to keep in mind regrading some of the dietary factors. For example calcium consumption was calculated based on consumption of dairy products, however calcium can be aquired from other sources including plant-based sources. However in these data, the influence of plant-based consumption of calcium was also not accounted for, nor was supplementation through vitamin sources.

What are the data?

We will be using data that we requested from the GBD about dietary intake, as well as the guideline data about optimal consumption amounts for different foods contained within the PDF of the article. We have two csv files. The first one includes consumption levels at the global level and for different countries for all ages combined.

Looking at the csv file in excel:

Here you can see that the data contains mean consumption values for both men and women in various countries at the national level in 2017 for various foods that may be problematic for health. The units for the food varies. So for example, the mean column in row that says “Diet low in fiber” indicates the average consumption level per person in that region and of that gender of fiber in grams per day.

The second csv file has similar data, but consumption levels for different age groups are separated.

The authors of this article obtained the data from a variety of sources including household budjet surveys and nutritional surveys regarding 24 hour recall of food consumption and 24 hour unrinary sodium analysis. The ata was derived from sales data from Euromonitor, data from the United Nations Food and Agriculture Organization (FAO), estimates about national availability of specific nutrients, from the Supply Utilazation Accounts(SUA), and the United States Department of Agriculture’s National Nutrition Database.

Note: While gender and sex are not actually binary, the data presented that is used in this analysis only contains data for groups of individuals described as men or women.

Data Import

Let’s import our data into R now so that we can explore the data further.

diet_data <-readr::read_csv(here("docs", "dietary_risk_exposure_all_ages_2017.csv"))
sep_age_diet_data <-read_csv(here("docs", "dietary_risk_exposure_sep_ages_2017.csv"))

First let’s just get a general sense of our data. We can do that using the glimpse() function of the dplyr package (it is also in the tibble package).

dplyr::glimpse(diet_data)

Observations: 5,880
Variables: 11
$ year_id        <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2…
$ location_name  <chr> "Global", "Global", "Global", "Global", "Global",…
$ rei_id         <dbl> 111, 111, 112, 112, 113, 113, 114, 114, 115, 115,…
$ rei_name       <chr> "Diet low in fruits", "Diet low in fruits", "Diet…
$ age_group_name <chr> "All Available Ages", "All Available Ages", "All …
$ sex            <chr> "Male", "Female", "Male", "Female", "Male", "Fema…
$ parameter      <chr> "continuous", "continuous", "continuous", "contin…
$ mean           <dbl> 8.703202e+01, 9.953182e+01, 1.975187e+02, 1.83233…
$ upper          <dbl> 8.988734e+01, 1.023450e+02, 2.083679e+02, 1.92110…
$ lower          <dbl> 8.452347e+01, 9.672140e+01, 1.887154e+02, 1.75401…
$ unit           <chr> "g/day", "g/day", "g/day", "g/day", "g/day", "g/d…

glimpse(sep_age_diet_data)

Observations: 88,200
Variables: 11
$ year_id        <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2…
$ location_name  <chr> "Global", "Global", "Global", "Global", "Global",…
$ rei_name       <chr> "Diet low in fruits", "Diet low in fruits", "Diet…
$ age_group_id   <dbl> 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 30, 3…
$ age_group_name <chr> "25 to 29", "30 to 34", "35 to 39", "40 to 44", "…
$ sex            <chr> "Male", "Male", "Male", "Male", "Male", "Male", "…
$ parameter      <chr> "continuous", "continuous", "continuous", "contin…
$ mean           <dbl> 68.54567, 72.62209, 76.73157, 81.03910, 88.79202,…
$ upper          <dbl> 75.16806, 80.00433, 84.73205, 89.52747, 100.94301…
$ lower          <dbl> 63.17081, 66.78018, 71.04691, 73.82734, 80.40971,…
$ unit           <chr> "g/day", "g/day", "g/day", "g/day", "g/day", "g/d…

Here we can tell that the sep_age_diet_data is much larger than the diet_data. There are 88,200 rows! The diet_data has only 5,880 rows. However, both files appear to have the same column structure with 11 variables each.

The skim() function of the skimr package is also really helpful for getting a general sense of your data.

skim(diet_data)

Data summary

Name	diet_data
Number of rows	5880
Number of columns	11
_______________________
Column type frequency:
character	6
numeric	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
location_name	0	1	4	32	0	196	0
rei_name	0	1	16	39	0	15	0
age_group_name	0	1	18	18	0	1	0
sex	0	1	4	6	0	2	0
parameter	0	1	10	10	0	1	0
unit	0	1	5	11	0	2	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
year_id	0	1	2017.00	0.00	2017	2017.00	2017.00	2017.00	2017.00	▁▁▇▁▁
rei_id	0	1	133.67	53.93	111	114.00	118.00	123.00	333.00	▇▁▁▁▁
mean	0	1	34.81	60.13	0	0.45	8.38	42.66	566.69	▇▁▁▁▁
upper	0	1	37.94	65.30	0	0.49	9.22	46.54	624.23	▇▁▁▁▁
lower	0	1	31.93	55.40	0	0.42	7.52	39.21	513.22	▇▁▁▁▁

Notice how there is a column about the values that are missing. It looks like our data is very complete and we do not have any missing data. We also get a sense about the size of our data.

The n_unqiue column shows us the number of unqiue values for each of our columns.

Let’s take a look at sep_age_diet_data.

skim(sep_age_diet_data)

Data summary

Name	sep_age_diet_data
Number of rows	88200
Number of columns	11
_______________________
Column type frequency:
character	6
numeric	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
location_name	0	1	4	32	0	196	0
rei_name	0	1	16	39	0	15	0
age_group_name	0	1	7	8	0	15	0
sex	0	1	4	6	0	2	0
parameter	0	1	10	10	0	1	0
unit	0	1	5	11	0	2	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
year_id	0	1	2017.00	0.00	2017	2017.00	2017.00	2017.00	2017.00	▁▁▇▁▁
age_group_id	0	1	32.87	54.46	10	13.00	17.00	30.00	235.00	▇▁▁▁▁
mean	0	1	39.84	63.86	0	0.56	10.80	52.63	604.01	▇▁▁▁▁
upper	0	1	48.07	78.38	0	0.70	13.55	59.96	806.98	▇▁▁▁▁
lower	0	1	33.46	54.54	0	0.42	8.46	44.87	639.06	▇▁▁▁▁

We can see that there are many more rows in this dataset.

Let’s take a look at the different dietary risk factors considered. To do this we will use the distinct() function of the dplyr package.

This function grabs only the distinct or unique rows from a given variable (rei_name, in our case) of a given data frame (diet_data, in our case).

#distinct(tibble_name, column_name)
  dplyr::distinct(diet_data, rei_name)

# A tibble: 15 x 1
   rei_name                               
   <chr>                                  
 1 Diet low in fruits                     
 2 Diet low in vegetables                 
 3 Diet low in whole grains               
 4 Diet low in nuts and seeds             
 5 Diet low in milk                       
 6 Diet high in red meat                  
 7 Diet high in processed meat            
 8 Diet high in sugar-sweetened beverages 
 9 Diet low in fiber                      
10 Diet low in seafood omega-3 fatty acids
11 Diet low in polyunsaturated fatty acids
12 Diet high in trans fatty acids         
13 Diet high in sodium                    
14 Diet low in calcium                    
15 Diet low in legumes                    

We will be using the %>% pipe for sequential steps in our code later on. This will make more sense when we have multiple sequential steps using the same data object.

We could do the same code as above using this notation. For example we first grab the diet_data, then we select the distinct values of the rei_name variable.

diet_data %>%
  distinct(rei_name)

# A tibble: 15 x 1
   rei_name                               
   <chr>                                  
 1 Diet low in fruits                     
 2 Diet low in vegetables                 
 3 Diet low in whole grains               
 4 Diet low in nuts and seeds             
 5 Diet low in milk                       
 6 Diet high in red meat                  
 7 Diet high in processed meat            
 8 Diet high in sugar-sweetened beverages 
 9 Diet low in fiber                      
10 Diet low in seafood omega-3 fatty acids
11 Diet low in polyunsaturated fatty acids
12 Diet high in trans fatty acids         
13 Diet high in sodium                    
14 Diet low in calcium                    
15 Diet low in legumes                    

Ok, so that gives us an idea of what dietary factors we can explore.

Let’s see if the location_name values are the same between both csv files. To do this we will use the setequal() function of dplyr.

dplyr::setequal(
  distinct(diet_data,location_name), 
  distinct(sep_age_diet_data, location_name)) 

[1] TRUE

Ok, we got the value of TRUE, so it looks like the same locations are in both files.

Note: In this case were comparing two different objects so using the pipe is not as useful.

Let’s take a look at the locations included in the data.

#scroll through the output!
sep_age_diet_data %>%
 distinct(location_name)%>%
  pull()

  [1] "Global"                           "China"                           
  [3] "North Korea"                      "Taiwan"                          
  [5] "Cambodia"                         "Indonesia"                       
  [7] "Laos"                             "Malaysia"                        
  [9] "Maldives"                         "Myanmar"                         
 [11] "Philippines"                      "Sri Lanka"                       
 [13] "Thailand"                         "Timor-Leste"                     
 [15] "Vietnam"                          "Fiji"                            
 [17] "Kiribati"                         "Marshall Islands"                
 [19] "Federated States of Micronesia"   "Papua New Guinea"                
 [21] "Samoa"                            "Solomon Islands"                 
 [23] "Tonga"                            "Vanuatu"                         
 [25] "Armenia"                          "Azerbaijan"                      
 [27] "Georgia"                          "Kazakhstan"                      
 [29] "Kyrgyzstan"                       "Mongolia"                        
 [31] "Tajikistan"                       "Turkmenistan"                    
 [33] "Uzbekistan"                       "Albania"                         
 [35] "Bosnia and Herzegovina"           "Bulgaria"                        
 [37] "Croatia"                          "Czech Republic"                  
 [39] "Hungary"                          "Macedonia"                       
 [41] "Montenegro"                       "Poland"                          
 [43] "Romania"                          "Serbia"                          
 [45] "Slovakia"                         "Slovenia"                        
 [47] "Belarus"                          "Estonia"                         
 [49] "Latvia"                           "Lithuania"                       
 [51] "Moldova"                          "Russian Federation"              
 [53] "Ukraine"                          "Brunei"                          
 [55] "Japan"                            "South Korea"                     
 [57] "Singapore"                        "Australia"                       
 [59] "New Zealand"                      "Andorra"                         
 [61] "Austria"                          "Belgium"                         
 [63] "Cyprus"                           "Denmark"                         
 [65] "Finland"                          "France"                          
 [67] "Germany"                          "Greece"                          
 [69] "Iceland"                          "Ireland"                         
 [71] "Israel"                           "Italy"                           
 [73] "Luxembourg"                       "Malta"                           
 [75] "Netherlands"                      "Norway"                          
 [77] "Portugal"                         "Spain"                           
 [79] "Sweden"                           "Switzerland"                     
 [81] "United Kingdom"                   "Argentina"                       
 [83] "Chile"                            "Uruguay"                         
 [85] "Canada"                           "United States"                   
 [87] "Antigua and Barbuda"              "The Bahamas"                     
 [89] "Barbados"                         "Belize"                          
 [91] "Cuba"                             "Dominica"                        
 [93] "Dominican Republic"               "Grenada"                         
 [95] "Guyana"                           "Haiti"                           
 [97] "Jamaica"                          "Saint Lucia"                     
 [99] "Saint Vincent and the Grenadines" "Suriname"                        
[101] "Trinidad and Tobago"              "Bolivia"                         
[103] "Ecuador"                          "Peru"                            
[105] "Colombia"                         "Costa Rica"                      
[107] "El Salvador"                      "Guatemala"                       
[109] "Honduras"                         "Mexico"                          
[111] "Nicaragua"                        "Panama"                          
[113] "Venezuela"                        "Brazil"                          
[115] "Paraguay"                         "Algeria"                         
[117] "Bahrain"                          "Egypt"                           
[119] "Iran"                             "Iraq"                            
[121] "Jordan"                           "Kuwait"                          
[123] "Lebanon"                          "Libya"                           
[125] "Morocco"                          "Palestine"                       
[127] "Oman"                             "Qatar"                           
[129] "Saudi Arabia"                     "Syria"                           
[131] "Tunisia"                          "Turkey"                          
[133] "United Arab Emirates"             "Yemen"                           
[135] "Afghanistan"                      "Bangladesh"                      
[137] "Bhutan"                           "India"                           
[139] "Nepal"                            "Pakistan"                        
[141] "Angola"                           "Central African Republic"        
[143] "Congo"                            "Democratic Republic of the Congo"
[145] "Equatorial Guinea"                "Gabon"                           
[147] "Burundi"                          "Comoros"                         
[149] "Djibouti"                         "Eritrea"                         
[151] "Ethiopia"                         "Kenya"                           
[153] "Madagascar"                       "Malawi"                          
[155] "Mauritius"                        "Mozambique"                      
[157] "Rwanda"                           "Seychelles"                      
[159] "Somalia"                          "Tanzania"                        
[161] "Uganda"                           "Zambia"                          
[163] "Botswana"                         "Lesotho"                         
[165] "Namibia"                          "South Africa"                    
[167] "Swaziland"                        "Zimbabwe"                        
[169] "Benin"                            "Burkina Faso"                    
[171] "Cameroon"                         "Cape Verde"                      
[173] "Chad"                             "Cote d'Ivoire"                   
[175] "The Gambia"                       "Ghana"                           
[177] "Guinea"                           "Guinea-Bissau"                   
[179] "Liberia"                          "Mali"                            
[181] "Mauritania"                       "Niger"                           
[183] "Nigeria"                          "Sao Tome and Principe"           
[185] "Senegal"                          "Sierra Leone"                    
[187] "Togo"                             "American Samoa"                  
[189] "Bermuda"                          "Greenland"                       
[191] "Guam"                             "Northern Mariana Islands"        
[193] "Puerto Rico"                      "Virgin Islands, U.S."            
[195] "South Sudan"                      "Sudan"                           

OK, so there are global values, as well as values for 185 countries.

Let’s take a look at the data when we order it by the mean consumption rate column:

diet_data %>%
  arrange(mean) %>%
  glimpse()

Observations: 5,880
Variables: 11
$ year_id        <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2…
$ location_name  <chr> "Lebanon", "Lebanon", "Italy", "Turkey", "Kazakhs…
$ rei_id         <dbl> 123, 123, 123, 123, 123, 123, 123, 123, 123, 123,…
$ rei_name       <chr> "Diet high in trans fatty acids", "Diet high in t…
$ age_group_name <chr> "All Available Ages", "All Available Ages", "All …
$ sex            <chr> "Male", "Female", "Male", "Male", "Male", "Male",…
$ parameter      <chr> "continuous", "continuous", "continuous", "contin…
$ mean           <dbl> 0.000727793, 0.000819837, 0.001357473, 0.00141107…
$ upper          <dbl> 0.000892791, 0.001007998, 0.001491509, 0.00172748…
$ lower          <dbl> 0.000585340, 0.000670438, 0.001233583, 0.00113497…
$ unit           <chr> "%energy/day", "%energy/day", "%energy/day", "%en…

Ok, so it looks like people in Lebanon dont eat very many trans fatty acids.

Let’s also figure out how many values there are in each age group of the data that is separated by age.

sep_age_diet_data %>%
dplyr::count(age_group_name)

# A tibble: 15 x 2
   age_group_name     n
   <chr>          <int>
 1 25 to 29        5880
 2 30 to 34        5880
 3 35 to 39        5880
 4 40 to 44        5880
 5 45 to 49        5880
 6 50 to 54        5880
 7 55 to 59        5880
 8 60 to 64        5880
 9 65 to 69        5880
10 70 to 74        5880
11 75 to 79        5880
12 80 to 84        5880
13 85 to 89        5880
14 90 to 94        5880
15 95 plus         5880

That’s a lot of values!

Let’s look a bit deeper to try to understand why. We can use the count function again but get the number of values for each category within sex, age_group_name and location_name of the data.

sep_age_diet_data%>%
  count(sex,age_group_name, location_name)

# A tibble: 5,880 x 4
   sex    age_group_name location_name           n
   <chr>  <chr>          <chr>               <int>
 1 Female 25 to 29       Afghanistan            15
 2 Female 25 to 29       Albania                15
 3 Female 25 to 29       Algeria                15
 4 Female 25 to 29       American Samoa         15
 5 Female 25 to 29       Andorra                15
 6 Female 25 to 29       Angola                 15
 7 Female 25 to 29       Antigua and Barbuda    15
 8 Female 25 to 29       Argentina              15
 9 Female 25 to 29       Armenia                15
10 Female 25 to 29       Australia              15
# … with 5,870 more rows

Ok, so it looks like there are probably the consumption values for each of the different dietary factors(which there were 15 different factors) for each age group, for each gender, and for each country.

We can confirm this by filtering the data to one of the age groups, for a single gender, and for a single location.

sep_age_diet_data %>%
filter( sex == "Female",
        age_group_name == "25 to 29",
        location_name == "Afghanistan")

# A tibble: 15 x 11
   year_id location_name rei_name age_group_id age_group_name sex  
     <dbl> <chr>         <chr>           <dbl> <chr>          <chr>
 1    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
 2    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
 3    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
 4    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
 5    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
 6    2017 Afghanistan   Diet hi…           10 25 to 29       Fema…
 7    2017 Afghanistan   Diet hi…           10 25 to 29       Fema…
 8    2017 Afghanistan   Diet hi…           10 25 to 29       Fema…
 9    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
10    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
11    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
12    2017 Afghanistan   Diet hi…           10 25 to 29       Fema…
13    2017 Afghanistan   Diet hi…           10 25 to 29       Fema…
14    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
15    2017 Afghanistan   Diet lo…           10 25 to 29       Fema…
# … with 5 more variables: parameter <chr>, mean <dbl>, upper <dbl>,
#   lower <dbl>, unit <chr>

Let’s also get the data from the PDF of the paper so that we can calculate consumption of these dietary factors as percentage of daily requirement, which would be more interpretable.

We are interested in this table on page 3:

First let’s import the PDF using the pdftools package.

paper<-pdftools::pdf_text(here("docs", "Afshin et al. 2019 - Health effects of dietary risks in 195 countries,  ... 17 - a systematic analysis for the Global Burden of Disease Study 2017.pdf"))

We can use the base summary() function to get a sense of what the data looks like. By base we mean that these functions are part of the base package and are loaded automatically.Thus library(base) is not required.

summary(paper)

   Length     Class      Mode 
       15 character character 

#This is equivalent to the following, but this is unecessary:
#base::summary(paper)

We can see that we have 15 different character strings. Each one contains the text on each of the 15 different pages of the PDF.

Data Wrangling

Again, the table we are interested in is on the third page, so let’s grab just that portion of the PDF.

#Here we will select the 3rd value in the paper object
table <- paper[3]

summary(table)

   Length     Class      Mode 
        1 character character 

glimpse(table, nchar.max = 800)

 chr "                                                                                                                                                                                                    Articles\nin systolic blood pressure, and then estimated the                                           Disease-specific deaths and disability-adjusted\nrelationship between change in systolic blood pressure                                       life-years\nand disease outcomes.14                                                                      Data on disease-specific deaths and disability-adjusted\n                                                                                             life-years (DALYs) by age, sex, country, and year were\nOptimal level of intake        "| __truncated__

Here we can see that the table object now contains the text from the 3rd page as a single large character string. However the text is difficult to read beacuse of the column structure in the pdf. Now let’s try to grab just the text in the table.

One way to approach this is to split the string by some pattern that we notice in the table.

Only the capitalized form of the word “Diet” apears to be within the table, and is not present in the preceding text (altough “diet” is). All the rows of interest of the table appear to start with the word “Diet”.

Let’s use the str_split() function of the stringr package to split the data within the object called tableby the word “Diet”. Only lines from page 3 that contain the word Diet will be selected (and not “diet” as this function is case-sensitive). Each section of the text that contians “Diet” will be split into individual pieces everytime the world “Diet” occurs and the word itself will be removed.

In this case we are also using the magrittr assignment pipe or double pipe that looks like this %<>%. This allows us use the table data as input to the later steps but also reassign the output to the same data object name.

table %<>%
  stringr::str_split(pattern = 'Diet')

Using the base::summary() and dplyr::glimpse() function we can see that we created a list of the rows in the table that contain the word “Diet”. We can see that we start with the row that contains “Diet low in fruits”.

table %>%
 summary()

     Length Class  Mode     
[1,] 17     -none- character

table %>%
  glimpse()

List of 1
 $ : chr [1:17] "                                                                                                               "| __truncated__ " low in fruits                  Mean daily consumption of fruits (fresh, frozen, cooked, canned, or dried fruit"| __truncated__ " low in vegetables              Mean daily consumption of vegetables (fresh, frozen, cooked, canned, or dried v"| __truncated__ " low in legumes                 Mean daily consumption of legumes (fresh, frozen, cooked, canned, or dried legu"| __truncated__ ...

RStudio creates really cheatsheets like this one which shows you all the major functions in stringr. You can download others here.

You can see that we could have also used the str_split_fixed() function which would also separate the substrings into different columns of a matrix, however we would need to know the number of substrings or pieces that we would like returned.

For more information about str_split() see here.

Let’s separte the values within the list using the base unlist function, this will allow us to easily select the different substrings within the object called table.

table %<>%
  unlist()

It’s important to realize that the first split will split the text before the first occurance of Data as the first value in the output. We could use the first() function of the dplyr package to look at this value. However, we will suppress the output as this is quite large.

dplyr::first(table)

Instead we can take a look at the second element of the list. using the nth() function of dplyr.

nth(table, 2)

[1] " low in fruits                  Mean daily consumption of fruits (fresh, frozen, cooked, canned, or dried fruits, excluding              250 g (200–300) per day                           94·9\n                                       fruit juices and salted or pickled fruits)\n   "

Indeed this looks like the first row of interest in our table:

Using the last() and the nth() functions of the dplyr package we can take a look at the last values of the list.

#to see the second to last value we can use nth()
#the -2 specifies that we want the second to last value
#-3 would be third to last and -1 would be the last value
dplyr::nth(table, -2)

[1] " high in sodium                 24 h urinary sodium measured in g per day                                                                   3 g (1–5) per day*                             26·2\n  *To reflect the uncertainty in existing evidence on optimal level of intake for sodium, 1–5 g per day was considered as the uncertainty range for the optimal level of sodium where less than 2·3 g per day is the\n  intake level of sodium associated with the lowest level of blood pressure in randomised controlled trials and 4–5 g per day is the level of sodium intake associated with the lowest risk of cardiovascular disease in\n  observational studies.\n  Table: "

#to see the very last value we can use last()
dplyr::last(table)

[1] "ary risk factor exposure definitions, optimal level, and data representativeness index, 1990–2017\nwww.thelancet.com Published online April 3, 2019 http://dx.doi.org/10.1016/S0140-6736(19)30041-8                                                                                                                         3\n"

Therefore, we dont need this part of the table or the text before the table if we just want the consumption reccomendations.

So we will select the 2nd through the second to last of the substrings. Since we have 17 substrings, we will select the 2nd through the 16th. However a better way to do this rather than selecting by index, would be to select phrases that are unique to the text within the table that we want. We will use the str_subset() function of stringr package to select the table rows with consumption guidlines. Most of the rows have the phrase" Mean daily consumption“, however, there are other phrases for some of the rows, including”Mean daily intake" and “24 h sodium”

# one could subset the table like this:
#table <- table[2:16]

table %<>%
str_subset(pattern = "Mean daily consumption|Mean daily intake|24 h")

Notice that we speparate the different patterns to look for using vertical bar character “|” and that all of the patterns are within quotation marks together.

Question opportunity:

1. What other string patterns could you use to subset the rows of the table that we want?

2. Why might it be better to subset based on the text rather than the index?

Now the first row is what we want:

first(table)

[1] " low in fruits                  Mean daily consumption of fruits (fresh, frozen, cooked, canned, or dried fruits, excluding              250 g (200–300) per day                           94·9\n                                       fruit juices and salted or pickled fruits)\n   "

And the last row is what we want:

last(table)

[1] " high in sodium                 24 h urinary sodium measured in g per day                                                                   3 g (1–5) per day*                             26·2\n  *To reflect the uncertainty in existing evidence on optimal level of intake for sodium, 1–5 g per day was considered as the uncertainty range for the optimal level of sodium where less than 2·3 g per day is the\n  intake level of sodium associated with the lowest level of blood pressure in randomised controlled trials and 4–5 g per day is the level of sodium intake associated with the lowest risk of cardiovascular disease in\n  observational studies.\n  Table: "

Notice that there the decimal points from the pdf are being recognized as an interpunct instead of a period or decimal. An interpunct is a centered dot, as opposed to a period or decimal that is aligned to the bottom of the line.

The interpunct was previously used to separate words in certain languages, like ancient Latin.

[source]

You can produce an interpunct on a mac like this:

[source]

It is important to replace these for later when we want these values to be converted from character strings to numeric. We will again use the stringr package. This time we will use the str_replace_all() function which replaces all instances of a pattern in an individual string. In this case we want to replace all instances of the interpunct with a decimal point.

table %<>%
  stringr::str_replace_all( pattern = "·", 
                            replacement = ".")

Now we will try to split the strings for each row based on the presence of 2 spaces to create the columns of the table, as there appears to be larger than a space between the columns to create substrings. The substrings will be separated by quotes.

The second page of the stringr cheetsheet has more information about using “Special Characters” in stringr. For example \\s is interpreted as a space as the \\ indicates that the s should be interpreted as a special character and not simply the letter s. The {2,} indicates 2 or more spaces, while {2} would indicate exactly 2 spaces.

table_split <- str_split(string=table, 
                         pattern= "\\s{2,}")
glimpse(table_split) #scroll the output!

List of 15
 $ : chr [1:6] " low in fruits" "Mean daily consumption of fruits (fresh, frozen, cooked, canned, or dried fruits, excluding" "250 g (200–300) per day" "94.9" ...
 $ : chr [1:7] " low in vegetables" "Mean daily consumption of vegetables (fresh, frozen, cooked, canned, or dried vegetables," "360 g (290–430) per day" "94.9" ...
 $ : chr [1:5] " low in legumes" "Mean daily consumption of legumes (fresh, frozen, cooked, canned, or dried legumes)" "60 g (50–70) per day" "94.9" ...
 $ : chr [1:7] " low in whole grains" "Mean daily consumption of whole grains (bran, germ, and endosperm in their natural" "125 g (100–150) per day" "94.9" ...
 $ : chr [1:5] " low in nuts and seeds" "Mean daily consumption of nut and seed foods" "21 g (16–25) per day" "94.9" ...
 $ : chr [1:6] " low in milk" "Mean daily consumption of milk including non-fat, low-fat, and full-fat milk, excluding soy" "435 g (350–520) per day" "94.9" ...
 $ : chr [1:6] " high in red meat" "Mean daily consumption of red meat (beef, pork, lamb, and goat, but excluding poultry, fish," "23 g (18–27) per day" "94.9" ...
 $ : chr [1:6] " high in processed meat" "Mean daily consumption of meat preserved by smoking, curing, salting, or addition of" "2 g (0–4) per day" "36.9" ...
 $ : chr [1:6] " high in sugar-sweetened Mean daily consumption of beverages with ≥50 kcal per 226.8 serving, including carbonated" "3 g (0–5) per day" "36.9" "beverages" ...
 $ : chr [1:6] " low in fibre" "Mean daily intake of fibre from all sources including fruits, vegetables, grains, legumes, and" "24 g (19–28) per day" "94.9" ...
 $ : chr [1:5] " low in calcium" "Mean daily intake of calcium from all sources, including milk, yogurt, and cheese" "1.25 g (1.00–1.50) per day" "94.9" ...
 $ : chr [1:5] " low in seafood omega-3 Mean daily intake of eicosapentaenoic acid and docosahexaenoic acid" "250 mg (200–300) per day" "94.9" "fatty acids" ...
 $ : chr [1:7] " low in polyunsaturated" "Mean daily intake of omega-6 fatty acids from all sources, mainly liquid vegetable oils," "11% (9–13) of total daily energy" "94.9" ...
 $ : chr [1:6] " high in trans fatty acids" "Mean daily intake of trans fat from all sources, mainly partially hydrogenated vegetable oils" "0.5% (0.0–1.0) of total daily energy" "36.9" ...
 $ : chr [1:8] " high in sodium" "24 h urinary sodium measured in g per day" "3 g (1–5) per day*" "26.2" ...

If we look closely, we can see that the sugar-sweetened beverage and the seafood category had only one space between the first and second columns - the columns about the dietary category and the one that describes in more detail what the consumption suggestion is about.

The values for these two columns appear to be together still in the same substring for these two categories. There are no quotation marks adjacent to the word "Mean".

Here you can see how the next substring should have started with the word "Mean" by the new inclusion of a quotation mark ".

We can add an extra space in front of the word "Mean" for these particular categories and then try splitting again.

Since we orginally split based on 2 or more spaces, we can just add a space in front of the word “Mean” for all the table strings and then try subsetting again.

table%>%
str_which(pattern = "seafood|sugar")

[1]  9 12

table[9] <-stringr::str_replace(pattern = "Mean", 
                                replacement = " Mean", table[9])
table[12] <-stringr::str_replace(pattern ="Mean", 
                                 replacement =" Mean", table[12])
table_split <- str_split(table,pattern= "\\s{2,}")

We could also hust add a space in front of all the values of Mean in the table since the split was peformed based on 2 or more spaces. Thus the other elements in table would also be split just as before despite the additional space.

table<-table %>%
  stringr::str_replace(pattern ="Mean", 
                       replacement = " Mean")
table_split <- str_split(table,pattern= "\\s{2,}")

#scroll the output!
glimpse(table_split) 

List of 15
 $ : chr [1:6] " low in fruits" "Mean daily consumption of fruits (fresh, frozen, cooked, canned, or dried fruits, excluding" "250 g (200–300) per day" "94.9" ...
 $ : chr [1:7] " low in vegetables" "Mean daily consumption of vegetables (fresh, frozen, cooked, canned, or dried vegetables," "360 g (290–430) per day" "94.9" ...
 $ : chr [1:5] " low in legumes" "Mean daily consumption of legumes (fresh, frozen, cooked, canned, or dried legumes)" "60 g (50–70) per day" "94.9" ...
 $ : chr [1:7] " low in whole grains" "Mean daily consumption of whole grains (bran, germ, and endosperm in their natural" "125 g (100–150) per day" "94.9" ...
 $ : chr [1:5] " low in nuts and seeds" "Mean daily consumption of nut and seed foods" "21 g (16–25) per day" "94.9" ...
 $ : chr [1:6] " low in milk" "Mean daily consumption of milk including non-fat, low-fat, and full-fat milk, excluding soy" "435 g (350–520) per day" "94.9" ...
 $ : chr [1:6] " high in red meat" "Mean daily consumption of red meat (beef, pork, lamb, and goat, but excluding poultry, fish," "23 g (18–27) per day" "94.9" ...
 $ : chr [1:6] " high in processed meat" "Mean daily consumption of meat preserved by smoking, curing, salting, or addition of" "2 g (0–4) per day" "36.9" ...
 $ : chr [1:7] " high in sugar-sweetened" "Mean daily consumption of beverages with ≥50 kcal per 226.8 serving, including carbonated" "3 g (0–5) per day" "36.9" ...
 $ : chr [1:6] " low in fibre" "Mean daily intake of fibre from all sources including fruits, vegetables, grains, legumes, and" "24 g (19–28) per day" "94.9" ...
 $ : chr [1:5] " low in calcium" "Mean daily intake of calcium from all sources, including milk, yogurt, and cheese" "1.25 g (1.00–1.50) per day" "94.9" ...
 $ : chr [1:6] " low in seafood omega-3" "Mean daily intake of eicosapentaenoic acid and docosahexaenoic acid" "250 mg (200–300) per day" "94.9" ...
 $ : chr [1:7] " low in polyunsaturated" "Mean daily intake of omega-6 fatty acids from all sources, mainly liquid vegetable oils," "11% (9–13) of total daily energy" "94.9" ...
 $ : chr [1:6] " high in trans fatty acids" "Mean daily intake of trans fat from all sources, mainly partially hydrogenated vegetable oils" "0.5% (0.0–1.0) of total daily energy" "36.9" ...
 $ : chr [1:8] " high in sodium" "24 h urinary sodium measured in g per day" "3 g (1–5) per day*" "26.2" ...

Looks better!

We want just the first (the food category) and third column (the optimal consumption amount suggested) for each row in the table.

We can use the map function of the purrr package to accomplish this.

The map function allows us to perform the same action multiple times across each element within an object.

This following will allow us to select the 1st or 3rd substring from each element of the table object.

category <-map(table_split,1)
amount <-map(table_split,3)
head(category)

[[1]]
[1] " low in fruits"

[[2]]
[1] " low in vegetables"

[[3]]
[1] " low in legumes"

[[4]]
[1] " low in whole grains"

[[5]]
[1] " low in nuts and seeds"

[[6]]
[1] " low in milk"

head(amount)

[[1]]
[1] "250 g (200–300) per day"

[[2]]
[1] "360 g (290–430) per day"

[[3]]
[1] "60 g (50–70) per day"

[[4]]
[1] "125 g (100–150) per day"

[[5]]
[1] "21 g (16–25) per day"

[[6]]
[1] "435 g (350–520) per day"

Now we will create a tibble using this data. However, currently both category and amount are of class list. To create a tibble we need to unlist the data to create vectors.

class(category)

[1] "list"

category %<>%unlist()
amount %<>%unlist()
class(category)

[1] "character"

category

 [1] " low in fruits"             " low in vegetables"        
 [3] " low in legumes"            " low in whole grains"      
 [5] " low in nuts and seeds"     " low in milk"              
 [7] " high in red meat"          " high in processed meat"   
 [9] " high in sugar-sweetened"   " low in fibre"             
[11] " low in calcium"            " low in seafood omega-3"   
[13] " low in polyunsaturated"    " high in trans fatty acids"
[15] " high in sodium"           

amount

 [1] "250 g (200–300) per day"             
 [2] "360 g (290–430) per day"             
 [3] "60 g (50–70) per day"                
 [4] "125 g (100–150) per day"             
 [5] "21 g (16–25) per day"                
 [6] "435 g (350–520) per day"             
 [7] "23 g (18–27) per day"                
 [8] "2 g (0–4) per day"                   
 [9] "3 g (0–5) per day"                   
[10] "24 g (19–28) per day"                
[11] "1.25 g (1.00–1.50) per day"          
[12] "250 mg (200–300) per day"            
[13] "11% (9–13) of total daily energy"    
[14] "0.5% (0.0–1.0) of total daily energy"
[15] "3 g (1–5) per day*"                  

We could have done all of this at once in one command like this:

category <-unlist(map(table_split,1))
amount <-unlist(map(table_split,3))

Now we will create a tibble, which is an important data frame structure in the tidyverse which allows us to use other packages in the tidyverse with our data.

We will name our tibble columns now as we create our tibble using the tibble() function of both the tidyr and the tibble packages, as names are required in tibbles.

guidelines <-tibble::tibble(category = category,
                              amount = amount)
guidelines

# A tibble: 15 x 2
   category                     amount                              
   <chr>                        <chr>                               
 1 " low in fruits"             250 g (200–300) per day             
 2 " low in vegetables"         360 g (290–430) per day             
 3 " low in legumes"            60 g (50–70) per day                
 4 " low in whole grains"       125 g (100–150) per day             
 5 " low in nuts and seeds"     21 g (16–25) per day                
 6 " low in milk"               435 g (350–520) per day             
 7 " high in red meat"          23 g (18–27) per day                
 8 " high in processed meat"    2 g (0–4) per day                   
 9 " high in sugar-sweetened"   3 g (0–5) per day                   
10 " low in fibre"              24 g (19–28) per day                
11 " low in calcium"            1.25 g (1.00–1.50) per day          
12 " low in seafood omega-3"    250 mg (200–300) per day            
13 " low in polyunsaturated"    11% (9–13) of total daily energy    
14 " high in trans fatty acids" 0.5% (0.0–1.0) of total daily energy
15 " high in sodium"            3 g (1–5) per day*                  

Looking pretty good!

However, we want to separate the different amounts within the amount column.

Recall what the orginal table looked like:

Separating values within a variable

We can use the tidyr::separate() function to separate the data within the amount column into three new columns based on the optimal level and the optimal range. We can separate the values based on the open parantheses "(" and the long dash "–" characters.

# The first column will be called optimal
# It will contain the 1st part of the amount column data before the 1st underscore"("
# The 2nd column will be called lower
# It will contain the data after the "("
# The 3rd column will be called upper 
# It will contain the 2nd part of the data based on the "–"

guidelines%<>% 
  tidyr::separate(amount, 
                  c("optimal", "lower", "upper"),
                  sep ="[[(|–]]") 
head(guidelines)

# A tibble: 6 x 4
  category                 optimal  lower upper       
  <chr>                    <chr>    <chr> <chr>       
1 " low in fruits"         "250 g " 200   300) per day
2 " low in vegetables"     "360 g " 290   430) per day
3 " low in legumes"        "60 g "  50    70) per day 
4 " low in whole grains"   "125 g " 100   150) per day
5 " low in nuts and seeds" "21 g "  16    25) per day 
6 " low in milk"           "435 g " 350   520) per day

Let’s Also create a new variable/column in our tibble that indicates the direction that can be harmful for each dietary factor.

guidelines%<>%
  separate(category, c("direction", "food"), sep = " in ")
guidelines

# A tibble: 15 x 5
   direction food              optimal   lower upper                     
   <chr>     <chr>             <chr>     <chr> <chr>                     
 1 " low"    fruits            "250 g "  200   300) per day              
 2 " low"    vegetables        "360 g "  290   430) per day              
 3 " low"    legumes           "60 g "   50    70) per day               
 4 " low"    whole grains      "125 g "  100   150) per day              
 5 " low"    nuts and seeds    "21 g "   16    25) per day               
 6 " low"    milk              "435 g "  350   520) per day              
 7 " high"   red meat          "23 g "   18    27) per day               
 8 " high"   processed meat    "2 g "    0     4) per day                
 9 " high"   sugar-sweetened   "3 g "    0     5) per day                
10 " low"    fibre             "24 g "   19    28) per day               
11 " low"    calcium           "1.25 g " 1.00  1.50) per day             
12 " low"    seafood omega-3   "250 mg " 200   300) per day              
13 " low"    polyunsaturated   "11% "    9     13) of total daily energy 
14 " high"   trans fatty acids "0.5% "   0.0   1.0) of total daily energy
15 " high"   sodium            "3 g "    1     5) per day*               

If we wanted to remove the direction variable we could use the purrr::modify_at() function:

guidelines %>% purrr::modify_at("direction",~NULL)

Data cleaning with regular expressions

Ok, looking better, but we still need a bit of cleaning to remove symbols and extra words from the columns. Some of the extra symbols include: "%", ")" and the "*".

The "*" and the ")" are what we call metacharacters or regular expressions. These are characters that have special meanings.

Now we need the "\\" to indicate that we want these characters to be matched exactly and not interpreted as the meaning of the symbol.

See here for more info about regular expressions in R.

Also here we have a bit of an example using the str_count() function of stringr, which counts the number of instences of a character string. In this case we will look for individual characters but you could also search for words or phrases.

regextest<-readr::read_file(here("docs", "regEx.txt"))
regextest

[1] "Testing for ts or\ttabs can be tricky.(yes, it really can!*)\n"

str_count(regextest,"t")#notice this doesn't include the t in the tab

[1] 5

str_count(regextest,"\t") #search for tab

[1] 1

str_count(regextest,"\\t")#search for tab

[1] 1

# this will not work because r thinks this is part of the code itself
#str_count(regextest, ")") 
# this will not work because r thinks this is part of the code itself
#str_count(regextest, "\)")
str_count(regextest, "\\)") #this works!

[1] 1

# this also does not work
#str_count(regextest, "*")
# nor does this
#str_count(regextest, "\*")
str_count(regextest, "\\*")#this works!

[1] 1

We also want to make a unit variable so that we can make sure that our units are consistent later.

guidelines %>%
pull(optimal) 

 [1] "250 g "  "360 g "  "60 g "   "125 g "  "21 g "   "435 g "  "23 g "  
 [8] "2 g "    "3 g "    "24 g "   "1.25 g " "250 mg " "11% "    "0.5% "  
[15] "3 g "   

Notice that the values that are percentages dont have spaces between the number and the unit. We can separate the optimal values by a space or a percent symbol "%" using "|" to indicate that we want to separate by either. In this case we will lose the “%” and will need to add it back to those values.

guidelines%<>%
  separate(optimal, into =c("optimal", "unit"), sep = " |%", remove = FALSE)
guidelines

# A tibble: 15 x 6
   direction food              lower optimal unit  upper                   
   <chr>     <chr>             <chr> <chr>   <chr> <chr>                   
 1 " low"    fruits            200   250     g     300) per day            
 2 " low"    vegetables        290   360     g     430) per day            
 3 " low"    legumes           50    60      g     70) per day             
 4 " low"    whole grains      100   125     g     150) per day            
 5 " low"    nuts and seeds    16    21      g     25) per day             
 6 " low"    milk              350   435     g     520) per day            
 7 " high"   red meat          18    23      g     27) per day             
 8 " high"   processed meat    0     2       g     4) per day              
 9 " high"   sugar-sweetened   0     3       g     5) per day              
10 " low"    fibre             19    24      g     28) per day             
11 " low"    calcium           1.00  1.25    g     1.50) per day           
12 " low"    seafood omega-3   200   250     mg    300) per day            
13 " low"    polyunsaturated   9     11      ""    13) of total daily ener…
14 " high"   trans fatty acids 0.0   0.5     ""    1.0) of total daily ene…
15 " high"   sodium            1     3       g     5) per day*             

Great, so to now we will add “%” to the unit variable for the low in polyunsaturated and high in trans fatty acids rows.

First we need to replace the empty values with NA using the na_if() function of the dplyr package.

guidelines %<>%
na_if("")
guidelines

# A tibble: 15 x 6
   direction food              lower optimal unit  upper                   
   <chr>     <chr>             <chr> <chr>   <chr> <chr>                   
 1 " low"    fruits            200   250     g     300) per day            
 2 " low"    vegetables        290   360     g     430) per day            
 3 " low"    legumes           50    60      g     70) per day             
 4 " low"    whole grains      100   125     g     150) per day            
 5 " low"    nuts and seeds    16    21      g     25) per day             
 6 " low"    milk              350   435     g     520) per day            
 7 " high"   red meat          18    23      g     27) per day             
 8 " high"   processed meat    0     2       g     4) per day              
 9 " high"   sugar-sweetened   0     3       g     5) per day              
10 " low"    fibre             19    24      g     28) per day             
11 " low"    calcium           1.00  1.25    g     1.50) per day           
12 " low"    seafood omega-3   200   250     mg    300) per day            
13 " low"    polyunsaturated   9     11      <NA>  13) of total daily ener…
14 " high"   trans fatty acids 0.0   0.5     <NA>  1.0) of total daily ene…
15 " high"   sodium            1     3       g     5) per day*             

Then to replace the NA values, we can use the replace_na() function in the tidyr package and the mutate() function of dplyr to specify which values to replace, in this case the NA values within the variable unit. Essentially this variable gets reasigned with the new values, as we mostly think of the mutate() function as creating new variables.

guidelines %<>% 
  dplyr::mutate(unit = replace_na(unit, "%"))
guidelines %>%
  filter(unit == "%")

# A tibble: 2 x 6
  direction food              lower optimal unit  upper                    
  <chr>     <chr>             <chr> <chr>   <chr> <chr>                    
1 " low"    polyunsaturated   9     11      %     13) of total daily energy
2 " high"   trans fatty acids 0.0   0.5     %     1.0) of total daily ener…

Let’s also move unit to be the last column. We can use the select() and everything() functions of the dplyr package to do this.

guidelines %<>%
  select(-unit,everything())

Here you can see Hadley Wickham’s (Chief Scientist at RStudio) explanation for this behavior of select():

https://github.com/tidyverse/dplyr/issues/2838#issuecomment-306062800

To remove all of the remaining extra characters and words we will again use the stringr package. This time we will use the str_remove_all() function to remove all instances of these characters.

guidelines <-as_tibble(
  map(
    guidelines, str_remove_all,
    pattern = "\\) per day|\\) of total daily energy|\\*"))

Nice! that’s pretty clean but we can do a bit more.

Data type conversion

One of the next things to notice about our data is the character classes of our variables.

Notice that the optimal amounts of consumption are currently of class character as indicated by the <chr> just below the column names / variable names of the guidelines tibble:

guidelines

# A tibble: 15 x 6
   direction food              lower optimal upper unit 
   <chr>     <chr>             <chr> <chr>   <chr> <chr>
 1 " low"    fruits            200   250     300   g    
 2 " low"    vegetables        290   360     430   g    
 3 " low"    legumes           50    60      70    g    
 4 " low"    whole grains      100   125     150   g    
 5 " low"    nuts and seeds    16    21      25    g    
 6 " low"    milk              350   435     520   g    
 7 " high"   red meat          18    23      27    g    
 8 " high"   processed meat    0     2       4     g    
 9 " high"   sugar-sweetened   0     3       5     g    
10 " low"    fibre             19    24      28    g    
11 " low"    calcium           1.00  1.25    1.50  g    
12 " low"    seafood omega-3   200   250     300   mg   
13 " low"    polyunsaturated   9     11      13    %    
14 " high"   trans fatty acids 0.0   0.5     1.0   %    
15 " high"   sodium            1     3       5     g    

To convert these values to numeric we can use the mutate_at() function of the dplyr package.

The mutate_at() function allows us to perform a function on specific columns/variables within a tibble. We need to indicate which variables that we would like to convert using vars(). In this case if we look at the beginning of the guidelines tibble, we can see that optimal, lower and upper should be converted. As these three columns are sequential, we can simply put a : between optimal and upper to indicate that we want all the variables in between these columns to be converted.

guidelines%<>%
  mutate_at(vars(lower:upper), as.numeric)
guidelines

# A tibble: 15 x 6
   direction food              lower optimal upper unit 
   <chr>     <chr>             <dbl>   <dbl> <dbl> <chr>
 1 " low"    fruits              200  250    300   g    
 2 " low"    vegetables          290  360    430   g    
 3 " low"    legumes              50   60     70   g    
 4 " low"    whole grains        100  125    150   g    
 5 " low"    nuts and seeds       16   21     25   g    
 6 " low"    milk                350  435    520   g    
 7 " high"   red meat             18   23     27   g    
 8 " high"   processed meat        0    2      4   g    
 9 " high"   sugar-sweetened       0    3      5   g    
10 " low"    fibre                19   24     28   g    
11 " low"    calcium               1    1.25   1.5 g    
12 " low"    seafood omega-3     200  250    300   mg   
13 " low"    polyunsaturated       9   11     13   %    
14 " high"   trans fatty acids     0    0.5    1   %    
15 " high"   sodium                1    3      5   g    

Great! Now these variables are of class <dbl> (stands for double) which indicates that they are numeric. Here is a link for more info on numeric classes in R.

If we had not replaced the "·" interpunct values to a period conversion from character to numeric will be problematic and will result in NA values.

Data value reassignments

We seem to have lost the word "beverages" from the "sugar-sweetened beverages" category, as well as "fatty acids" from the "seafood omega 3 fatty acids", and the "polyunsaturated fatty acids" categories as the full category name was listed on two lines within the table. We would like to replace these values with the full name.

To select the food column we will show you several options. Only a couple will work well with reassigning the data in that particular variable within guidelines without assigning an intermediate data object. We will look using mutate_at(), pull(), select(), and brackets [,c("variable name")].

The bracket option and the select() option will grab a tibble (data frame) version of the food column out of guidelines. However we can’t start commands with select for assignments.

guidelines[,c("food")] #same output as select

# A tibble: 15 x 1
   food             
   <chr>            
 1 fruits           
 2 vegetables       
 3 legumes          
 4 whole grains     
 5 nuts and seeds   
 6 milk             
 7 red meat         
 8 processed meat   
 9 sugar-sweetened  
10 fibre            
11 calcium          
12 seafood omega-3  
13 polyunsaturated  
14 trans fatty acids
15 sodium           

select(guidelines, "food") # same output as brackets

# A tibble: 15 x 1
   food             
   <chr>            
 1 fruits           
 2 vegetables       
 3 legumes          
 4 whole grains     
 5 nuts and seeds   
 6 milk             
 7 red meat         
 8 processed meat   
 9 sugar-sweetened  
10 fibre            
11 calcium          
12 seafood omega-3  
13 polyunsaturated  
14 trans fatty acids
15 sodium           

pull() in contrast, will grab the vector version of the food data:

pull(guidelines, "food") # get character vector not a tibble

 [1] "fruits"            "vegetables"        "legumes"          
 [4] "whole grains"      "nuts and seeds"    "milk"             
 [7] "red meat"          "processed meat"    "sugar-sweetened"  
[10] "fibre"             "calcium"           "seafood omega-3"  
[13] "polyunsaturated"   "trans fatty acids" "sodium"           

The pull function can be very useful when combined with other functions (for example you typically want to use a vector with the str_replace() function), but just like select, we can’t start assignments with pull().

This is not possible and will result in an error:

select(guidelines, food) <- 
   str_replace( 
   pull(guidelines,"food"), 
   pattern = "sugar-sweetened", 
   replacement = "sugar-sweetened beverages")

This will only print the result, but not reassign the food variable values:

guidelines %>%
   pull(food)%>%
   str_replace( 
   pattern = "sugar-sweetened", 
   replacement = "sugar-sweetened beverages")

 [1] "fruits"                    "vegetables"               
 [3] "legumes"                   "whole grains"             
 [5] "nuts and seeds"            "milk"                     
 [7] "red meat"                  "processed meat"           
 [9] "sugar-sweetened beverages" "fibre"                    
[11] "calcium"                   "seafood omega-3"          
[13] "polyunsaturated"           "trans fatty acids"        
[15] "sodium"                   

Using select() would work as well to print the result (although the result structure is different):

guidelines %>%
   select(food)%>%
   str_replace( 
   pattern = "sugar-sweetened", 
   replacement = "sugar-sweetened beverages")

[1] "c(\"fruits\", \"vegetables\", \"legumes\", \"whole grains\", \"nuts and seeds\", \"milk\", \"red meat\", \"processed meat\", \"sugar-sweetened beverages\", \"fibre\", \"calcium\", \"seafood omega-3\", \"polyunsaturated\", \"trans fatty acids\", \"sodium\")"

Question opportunity:

Why do these commands not reassign the food variable values?

The bracket option is great alternative and allows us to reassign the values within guidelines easily

#Replacing "sugar-sweetened" with "sugar-sweetened beverages"
guidelines[,c("food")] <- 
  str_replace( 
  pull(guidelines,"food"), 
  pattern = "sugar-sweetened", 
  replacement = "sugar-sweetened beverages")

#Replacing "seafood omega-3" with"seafood omega-3 fatty acids"
guidelines[,c("food")] <- 
  str_replace( 
  pull(guidelines,"food"), 
  pattern = "seafood omega-3", 
  replacement = "seafood omega-3 fatty acids")

guidelines

# A tibble: 15 x 6
   direction food                        lower optimal upper unit 
   <chr>     <chr>                       <dbl>   <dbl> <dbl> <chr>
 1 " low"    fruits                        200  250    300   g    
 2 " low"    vegetables                    290  360    430   g    
 3 " low"    legumes                        50   60     70   g    
 4 " low"    whole grains                  100  125    150   g    
 5 " low"    nuts and seeds                 16   21     25   g    
 6 " low"    milk                          350  435    520   g    
 7 " high"   red meat                       18   23     27   g    
 8 " high"   processed meat                  0    2      4   g    
 9 " high"   sugar-sweetened beverages       0    3      5   g    
10 " low"    fibre                          19   24     28   g    
11 " low"    calcium                         1    1.25   1.5 g    
12 " low"    seafood omega-3 fatty acids   200  250    300   mg   
13 " low"    polyunsaturated                 9   11     13   %    
14 " high"   trans fatty acids               0    0.5    1   %    
15 " high"   sodium                          1    3      5   g    

Finally, the best option is probably the mutate_at() function from dplyr. In this case we need to include ~ in front of the function that we would like to use on the values in our food variables. We also include . as a replacement to reference the data that we want to use within str_replace() (which in this case is the food variable values of guidelines).

Notice we didn’t need this when we previously use mutate_at() with the as.numeric() function. This is becuase the str_replace() function requires us to specify what data we are using as one of the arguments, while as.numeric() does not.

#Replacing "polyunsaturated" with"polyunsaturated fatty acids"
guidelines%<>%
  mutate_at(vars(food),
  ~str_replace( 
  string = ., 
  pattern = "polyunsaturated", 
  replacement = "polyunsaturated fatty acids"))

guidelines

# A tibble: 15 x 6
   direction food                        lower optimal upper unit 
   <chr>     <chr>                       <dbl>   <dbl> <dbl> <chr>
 1 " low"    fruits                        200  250    300   g    
 2 " low"    vegetables                    290  360    430   g    
 3 " low"    legumes                        50   60     70   g    
 4 " low"    whole grains                  100  125    150   g    
 5 " low"    nuts and seeds                 16   21     25   g    
 6 " low"    milk                          350  435    520   g    
 7 " high"   red meat                       18   23     27   g    
 8 " high"   processed meat                  0    2      4   g    
 9 " high"   sugar-sweetened beverages       0    3      5   g    
10 " low"    fibre                          19   24     28   g    
11 " low"    calcium                         1    1.25   1.5 g    
12 " low"    seafood omega-3 fatty acids   200  250    300   mg   
13 " low"    polyunsaturated fatty acids     9   11     13   %    
14 " high"   trans fatty acids               0    0.5    1   %    
15 " high"   sodium                          1    3      5   g    

This might be considered a better option because it is more readible as to where the food data came from that we are replacing values within.

There is one last minor detail… the direction variable has leading spaces still. We can use str_trim() to fix that!

guidelines%<>%
  mutate_at(vars(direction), str_trim)

guidelines

# A tibble: 15 x 6
   direction food                        lower optimal upper unit 
   <chr>     <chr>                       <dbl>   <dbl> <dbl> <chr>
 1 low       fruits                        200  250    300   g    
 2 low       vegetables                    290  360    430   g    
 3 low       legumes                        50   60     70   g    
 4 low       whole grains                  100  125    150   g    
 5 low       nuts and seeds                 16   21     25   g    
 6 low       milk                          350  435    520   g    
 7 high      red meat                       18   23     27   g    
 8 high      processed meat                  0    2      4   g    
 9 high      sugar-sweetened beverages       0    3      5   g    
10 low       fibre                          19   24     28   g    
11 low       calcium                         1    1.25   1.5 g    
12 low       seafood omega-3 fatty acids   200  250    300   mg   
13 low       polyunsaturated fatty acids     9   11     13   %    
14 high      trans fatty acids               0    0.5    1   %    
15 high      sodium                          1    3      5   g    

OK! Now we know how much of each dietary factor we generally need for optimal health according to the guidelines used in this article.

We would like to see how the mean consumption rates for the different groups of people compared to the optimal intake guidelines.

One way we could do this is to calculate a consumption percentage of the optimal value.

To calculate this it would be helpful to put the guideline amounts with the average consumption rates into the same tibble, especially because the observed consumption data (diet_data and sep_age_diet_data) are very different dimensions from the guidelines data.

In order to create a tibble with our observed consumption rates with the suggested consumption rates, we will join our data using dplyr. In order to do so it is important that our different datasets have the same values. So let’s first assess if that is the case.

Comparing data

distinct(diet_data, rei_name)

# A tibble: 15 x 1
   rei_name                               
   <chr>                                  
 1 Diet low in fruits                     
 2 Diet low in vegetables                 
 3 Diet low in whole grains               
 4 Diet low in nuts and seeds             
 5 Diet low in milk                       
 6 Diet high in red meat                  
 7 Diet high in processed meat            
 8 Diet high in sugar-sweetened beverages 
 9 Diet low in fiber                      
10 Diet low in seafood omega-3 fatty acids
11 Diet low in polyunsaturated fatty acids
12 Diet high in trans fatty acids         
13 Diet high in sodium                    
14 Diet low in calcium                    
15 Diet low in legumes                    

select(guidelines, food)

# A tibble: 15 x 1
   food                       
   <chr>                      
 1 fruits                     
 2 vegetables                 
 3 legumes                    
 4 whole grains               
 5 nuts and seeds             
 6 milk                       
 7 red meat                   
 8 processed meat             
 9 sugar-sweetened beverages  
10 fibre                      
11 calcium                    
12 seafood omega-3 fatty acids
13 polyunsaturated fatty acids
14 trans fatty acids          
15 sodium                     

We can see that we need to remove the "Diet low in" and "Diet high in" phrases from the observed consumption data.

diet_data$rei_name<- diet_data$rei_name %>%
  str_remove( pattern = "Diet low in |Diet high in ")

sep_age_diet_data$rei_name<-sep_age_diet_data$rei_name %>%
  str_remove( pattern = "Diet low in |Diet high in ")

Also let’s double check that the two observed files have the same exact values for dietary factor names.

We can use the setequal() function from dplyr to check that the unique values for rei_name are the same for both diet_data and sep_age_diet_data.

setequal(distinct(diet_data, rei_name), 
         distinct(sep_age_diet_data, rei_name))

[1] TRUE

Great!

Note that the default of the set_equal function ignores the order of values in rows. So we still dont know if the order is the same.

We can check using the all_equal function of dplyr which reports back clues about what might be different if anything. Importantly we are including ignore_row_order = FALSE as the default is TRUE.

all_equal(distinct(diet_data, rei_name), 
          distinct(sep_age_diet_data, rei_name), 
          ignore_row_order = FALSE)

[1] TRUE

Looks like they are also in the same order.

Note that if any of the values are different all_equal() will first report this and will not report that the rows are in a different order.

Here is a toy example about how the three comparison functions(setequal(), all_equal() (also all.equal() for tbl_df), and setdiff()) work in dplyr.

It’s important to realize that row order is ignored by bothsetequal() and setdiff().

Now let’s compare two tibbles that have different row orders and differnt values.

Here are our tibbles to compare:

X <-tibble(test =c("A", "B", "AC", "D"))
Y <-tibble(test =c("A", "D", "AG", "B"))
X

# A tibble: 4 x 1
  test 
  <chr>
1 A    
2 B    
3 AC   
4 D    

Y

# A tibble: 4 x 1
  test 
  <chr>
1 A    
2 D    
3 AG   
4 B    

class(Y)

[1] "tbl_df"     "tbl"        "data.frame"

Since we are using tibbles, which are of class tbl_df we can use either all_equal or all.equal().

all_equal(X, Y, ignore_row_order = TRUE)

[1] "Rows in x but not y: 3. Rows in y but not x: 3. "

all_equal(X, Y, ignore_row_order = FALSE)

[1] "Rows in x but not y: 3. Rows in y but not x: 3. "

# doesnt report rows being different order
all.equal(X, Y, ignore_row_order = TRUE)

[1] "Rows in x but not y: 3. Rows in y but not x: 3. "

all.equal(X, Y, ignore_row_order = FALSE)

[1] "Rows in x but not y: 3. Rows in y but not x: 3. "

# doesnt report rows being different order

# Reports false indicating at least one difference
# Does not provide clues about what is different
setequal(X, Y)

[1] FALSE

#setdiff() tells us what is different
#setdiff() is dependent on the order of the objects compared
#This reports what is unique to X
setdiff(X, Y) 

# A tibble: 1 x 1
  test 
  <chr>
1 AC   

#This reports what is unique to Y
setdiff(Y, X) 

# A tibble: 1 x 1
  test 
  <chr>
1 AG   

Now let’s make it so that only the order is different:

Y <-tibble(test =c("A", "D", "AC", "B"))
X

# A tibble: 4 x 1
  test 
  <chr>
1 A    
2 B    
3 AC   
4 D    

Y

# A tibble: 4 x 1
  test 
  <chr>
1 A    
2 D    
3 AC   
4 B    

all_equal(X, Y, ignore_row_order = TRUE)

[1] TRUE

all_equal(X, Y, ignore_row_order = FALSE)# reports diff order

[1] "Same row values, but different order"

# Remember setequal() ignores order!
# It reports no difference!
setequal(X, Y)

[1] TRUE

# Set diff also ignores order!
setdiff(X, Y) 

# A tibble: 0 x 1
# … with 1 variable: test <chr>

If we have different column/variable names this makes comparisons more challenging:

X <-tibble(test =c("A", "B", "AC", "D"))
Y <-tibble(test2 =c("A", "D", "AG", "B"))
# just reports that col names are diff
all_equal(X, Y, ignore_row_order = TRUE)

[1] "Cols in y but not x: `test2`. " "Cols in x but not y: `test`. " 

# just reports that col names are diff
all_equal(X, Y, ignore_row_order = FALSE)

[1] "Cols in y but not x: `test2`. " "Cols in x but not y: `test`. " 

setequal(X, Y)#This works!

[1] FALSE

#setdiff(X, Y) # This will not work

Ok, let’s keep going with our data.

How similar are the guidelines tibble and the observed consumption tibbles?

setequal(distinct(diet_data, rei_name), 
          select(guidelines, food))

[1] FALSE

Ok, looks like we have some different values.

Let’s use the setdiff function to get more information about what is different between the values.

setdiff(distinct(diet_data, rei_name), 
          select(guidelines, food))

:( That wont work. This is because setdiff() requires that the colnames are the same in the objects that we are comparing.

We can use the rename() of dplyr function to do this. We list the value that we want to change to first. We find “food” more intuitive so we are going to change “rei_name” to “food” for the diet_data and the sep_age_diet_data

diet_data %<>%
  dplyr::rename(food = rei_name)
sep_age_diet_data %<>%
  dplyr::rename(food = rei_name)

setdiff(distinct(diet_data, food), 
          select(guidelines, food), 
          ignore_row_order = TRUE)

Great, now we know that the fiber value appears to be different between the two.

setdiff(select(guidelines, food),
        distinct(diet_data, food))

# A tibble: 1 x 1
  food 
  <chr>
1 fibre

Changing the order of the objects, we can see that in the table from the article that we used to create guidelines, “fibre” the British spelling is used in contrast to the dataset which uses the American spelling “fiber”.

Let’s stick with the American spelling, so we will replace "fibre" in the guideline tibble. Again, we have two options for doing this:

# guidelines[,c("food")] <- 
#   str_replace( 
#   pull(guidelines,"food"), 
#   pattern = "fibre", 
#   replacement = "fiber")

guidelines%<>%
  mutate_at(vars(food),
  ~str_replace( 
  string = ., 
  pattern = "fibre", 
  replacement = "fiber"))

guidelines %>%
  filter(food == "fiber")

# A tibble: 1 x 6
  direction food  lower optimal upper unit 
  <chr>     <chr> <dbl>   <dbl> <dbl> <chr>
1 low       fiber    19      24    28 g    

Now let’s check again to see that our food values match between the guidelines and the observed consumption data tibbles.

setdiff(select(guidelines, food),
        distinct(diet_data, food))

# A tibble: 0 x 1
# … with 1 variable: food <chr>

setdiff(select(guidelines, food),
        distinct(sep_age_diet_data, food))

# A tibble: 0 x 1
# … with 1 variable: food <chr>

Great! There are no differences :)

Joining data

Now we can put our guideline data together with the diet_data and the sep_age_diet_data.

Remeber that the food data in our guidelines tibble is not necessarily in the same order as that of the consumption data tibbles. Thus this could be a problem if we decided to expand the guidelines rows (to repeat for the number of fruit observations etc.) and add them to our observed consumption tibbles.

In that case we could use the arrange() function of dplyr to sort the data alphabetically.

However, we will instead use a joining function of dplyr. These functions combine the data together based on common values. There are a variety of options.

In our case we would like to retain all of the values of diet_data and sep_age_diet_data. We would like to add new columns of values to these tibbles that correspond to the guideline information about amounts of consumption for each food type in the guidelines tibble. We shouldn’t have any values of food in guidelines that dont match, so we will not get any NA values. Therefore, in our case any of the mutating join functions should result in the same output.

It’s important to check if we have any overlapping variable names before we join the data. We can use the base R function names() and the intersect() function of the dplyr package for this.

dplyr::intersect(names(diet_data), 
          names(guidelines))

[1] "food"  "upper" "lower" "unit" 

So it looks like the "upper" , "lower" and "unit" variable names are overlapping. Therefore, to distinguish the names later we will rename the guideline "upper" , "lower" and "unit" variables.

We will again use the rename function from the dplyr package. We can list multiple variables to rename and separate each with a comma.

guidelines %<>%
  rename(upper_optimal = upper, 
         lower_optimal = lower,
          unit_optimal = unit)

guidelines

# A tibble: 15 x 6
   direction food          lower_optimal optimal upper_optimal unit_optimal
   <chr>     <chr>                 <dbl>   <dbl>         <dbl> <chr>       
 1 low       fruits                  200  250            300   g           
 2 low       vegetables              290  360            430   g           
 3 low       legumes                  50   60             70   g           
 4 low       whole grains            100  125            150   g           
 5 low       nuts and see…            16   21             25   g           
 6 low       milk                    350  435            520   g           
 7 high      red meat                 18   23             27   g           
 8 high      processed me…             0    2              4   g           
 9 high      sugar-sweete…             0    3              5   g           
10 low       fiber                    19   24             28   g           
11 low       calcium                   1    1.25           1.5 g           
12 low       seafood omeg…           200  250            300   mg          
13 low       polyunsatura…             9   11             13   %           
14 high      trans fatty …             0    0.5            1   %           
15 high      sodium                    1    3              5   g           

It’s also a good idea to check our units to make sure they are the same for both guidelines and the observed consumption tibbles: diet_and_guidelines and all_age_diet_and_guidelines.

Let’s take a look with the count() function of dplyr.

bind_cols(dietdata = count(diet_data, unit, food), 
       sepagedietdata = count(sep_age_diet_data,unit,food),
       guideline= count(guidelines, unit_optimal, food))

# A tibble: 15 x 9
   unit   food          n unit1  food1       n1 unit_optimal food2       n2
   <chr>  <chr>     <int> <chr>  <chr>    <int> <chr>        <chr>    <int>
 1 %ener… polyunsa…   392 %ener… polyuns…  5880 %            polyuns…     1
 2 %ener… trans fa…   392 %ener… trans f…  5880 %            trans f…     1
 3 g/day  calcium     392 g/day  calcium   5880 g            calcium      1
 4 g/day  fiber       392 g/day  fiber     5880 g            fiber        1
 5 g/day  fruits      392 g/day  fruits    5880 g            fruits       1
 6 g/day  legumes     392 g/day  legumes   5880 g            legumes      1
 7 g/day  milk        392 g/day  milk      5880 g            milk         1
 8 g/day  nuts and…   392 g/day  nuts an…  5880 g            nuts an…     1
 9 g/day  processe…   392 g/day  process…  5880 g            process…     1
10 g/day  red meat    392 g/day  red meat  5880 g            red meat     1
11 g/day  seafood …   392 g/day  seafood…  5880 g            sodium       1
12 g/day  sodium      392 g/day  sodium    5880 g            sugar-s…     1
13 g/day  sugar-sw…   392 g/day  sugar-s…  5880 g            vegetab…     1
14 g/day  vegetabl…   392 g/day  vegetab…  5880 g            whole g…     1
15 g/day  whole gr…   392 g/day  whole g…  5880 mg           seafood…     1

We can see that the only potential issue is the seafood omega-3 fatty acids data which is in g/day for the observed data(diet_data and all_age_diet_and_guidelines), but the unit is mg/day in the guidelines data.

We can account for this by dividing the guideline seafood omega-3 fatty acids data by 1000 to convert it to grams from milligrams.

To do this we will use the if_else() function in the dplyr package. This allows us to specify a condition (in this case if the unit is "mg"), as well as values if this is codition is met (true), or if the condition is not met (false). In the following we mutate the values in each of the guideline numeric columns (lower, optimal and upper) one at a time. When we refer to lower for example we refer to the values in the column/varaible. So if the condition is not met, then the original value is retained. We will also replace "mg" with "g" after everything is converted to grams.

guidelines%<>% mutate(lower_optimal = if_else(
  condition = unit_optimal=="mg", 
  true = lower_optimal/1000, 
  false = lower_optimal))

guidelines%<>% mutate(optimal = if_else(
  condition = unit_optimal=="mg", 
  true = optimal/1000, 
  false = optimal))

guidelines%<>% mutate(upper_optimal = if_else(
  condition = unit_optimal=="mg", 
  true = upper_optimal/1000, 
  false = upper_optimal))


guidelines%<>% mutate(unit_optimal = if_else(
  condition = unit_optimal=="mg", 
  true = "g", 
  false = unit_optimal))


#or this:
# guidelines%<>%
#   mutate_at(vars(unit_optimal),
#   ~str_replace( 
#   string = ., 
#   pattern = "mg", 
#   replacement = "g"))


guidelines

# A tibble: 15 x 6
   direction food          lower_optimal optimal upper_optimal unit_optimal
   <chr>     <chr>                 <dbl>   <dbl>         <dbl> <chr>       
 1 low       fruits                200    250            300   g           
 2 low       vegetables            290    360            430   g           
 3 low       legumes                50     60             70   g           
 4 low       whole grains          100    125            150   g           
 5 low       nuts and see…          16     21             25   g           
 6 low       milk                  350    435            520   g           
 7 high      red meat               18     23             27   g           
 8 high      processed me…           0      2              4   g           
 9 high      sugar-sweete…           0      3              5   g           
10 low       fiber                  19     24             28   g           
11 low       calcium                 1      1.25           1.5 g           
12 low       seafood omeg…           0.2    0.25           0.3 g           
13 low       polyunsatura…           9     11             13   %           
14 high      trans fatty …           0      0.5            1   %           
15 high      sodium                  1      3              5   g           

THIS IS WORK to try to mutate all at once… can delete or keep working on…

#guidelines %>% mutate_if(str_detect(guidelines$unit, "mg") & is.numeric, ~.+2)
#https://stackoverflow.com/questions/51877611/can-i-combine-a-dplyr-mutate-at-mutate-if-statement
#https://github.com/tidyverse/dplyr/issues/631
#WORK
#if numeric and unit = mg??
# guidelines %>% 
#    mutate_at(filter(guidelines,food == "seafood omega-3 fatty acids"),vars(lower:upper), funs(./1000))
# 
# guidelines %>% 
#    if_else(unit == "mg", mutate_at(list(lower:upper),if_else(unit == "mg", true = funs(./1000), false = .)))
 # mutate_at(filter(guidelines,food == "seafood omega-3 fatty #acids"),vars(lower:upper), funs(./1000))
# guidelines %>% 
#    filter(food == "seafood omega-3 fatty acids")%>%
#    mutate_at(vars(lower:upper), funs(./1000))
# 
# guidelines %>% 
#    if_else(unit == "mg", mutate_at(list(lower:upper),if_else(unit == "mg", true = funs(./1000), false = .)))
# guidelines%>%
#   mutate_at(vars(lower:upper), funs(ifelse(unit=="mg", TRUE/1000, TRUE)))
# 
# 
# mutate(cars, dist = ifelse(speed==4, dist*100, dist))
# 
# 
# guidelines[c("seafood omega-3 fatty acids"),] 
# guidelines[c("seafood omega-3 fatty acids"),]  <-guidelines %>% 
#   filter(food == "seafood omega-3 fatty acids")%>%
#   mutate_at(vars(lower:upper), funs(10*.))
# 
# cars<-mutate(cars, dist2 = dist)
# mutate(cars, dist = ifelse(speed==4, dist*100, dist))

guidelines %>% 
    filter(food == "seafood omega-3 fatty acids")%>%
    mutate_at(vars(lower_optimal:upper_optimal), funs(./1000))

# A tibble: 1 x 6
  direction food           lower_optimal optimal upper_optimal unit_optimal
  <chr>     <chr>                  <dbl>   <dbl>         <dbl> <chr>       
1 low       seafood omega…        0.0002 0.00025      0.000300 g           

# type <- c(1:4)
# year1 <- c(1:4)
# year2 <- c(1:4)
# year3 <- c(1:4)
# data <- data.frame(type, year1, year2, year3)
# 
# 
# newdata <- data %>%
#   gather(., year, value, year1:year3) %>%
#   mutate(newvalue = ifelse(type > 2, value * 2, value)) %>%
#   select(-value) %>%
#   spread(., year, newvalue)
# 
# Df %>% 
#   mutate_all(
#     funs(case_when(
#     . == "1"  ~ 1*.,
#     . == "2"  ~ 2*.,
#     . == "3"  ~ 3*.))) %>%
# mutate(rowmax = pmax(!!!rlang::syms(names(.))))
# 
# df %>% mutate_at(vars(a:c), funs(ifelse(is.na(d) | d == 0, NA, .)))
# 
# 
# 
# msleep %>%
#   select(name, sleep_total) %>%
#   mutate(sleep_total_min = sleep_total * 60)

Now we are ready to join the data!

Agian, we would like to add new columns of values to diet_data and all_age_diet_and_guidelines that correspond to the guideline information about amounts of consumption for each food type in the guidelines tibble. So we will join the data based on the food variable values.

diet_and_guidelines <-diet_data %>%
  full_join(guidelines, by = "food" )

all_age_diet_and_guidelines <-sep_age_diet_data %>%
  full_join(guidelines, by = "food" )

glimpse(diet_and_guidelines)

Observations: 5,880
Variables: 16
$ year_id        <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2…
$ location_name  <chr> "Global", "Global", "Global", "Global", "Global",…
$ rei_id         <dbl> 111, 111, 112, 112, 113, 113, 114, 114, 115, 115,…
$ food           <chr> "fruits", "fruits", "vegetables", "vegetables", "…
$ age_group_name <chr> "All Available Ages", "All Available Ages", "All …
$ sex            <chr> "Male", "Female", "Male", "Female", "Male", "Fema…
$ parameter      <chr> "continuous", "continuous", "continuous", "contin…
$ mean           <dbl> 8.703202e+01, 9.953182e+01, 1.975187e+02, 1.83233…
$ upper          <dbl> 8.988734e+01, 1.023450e+02, 2.083679e+02, 1.92110…
$ lower          <dbl> 8.452347e+01, 9.672140e+01, 1.887154e+02, 1.75401…
$ unit           <chr> "g/day", "g/day", "g/day", "g/day", "g/day", "g/d…
$ direction      <chr> "low", "low", "low", "low", "low", "low", "low", …
$ lower_optimal  <dbl> 200.0, 200.0, 290.0, 290.0, 100.0, 100.0, 16.0, 1…
$ optimal        <dbl> 250.00, 250.00, 360.00, 360.00, 125.00, 125.00, 2…
$ upper_optimal  <dbl> 300.0, 300.0, 430.0, 430.0, 150.0, 150.0, 25.0, 2…
$ unit_optimal   <chr> "g", "g", "g", "g", "g", "g", "g", "g", "g", "g",…

glimpse(all_age_diet_and_guidelines)

Observations: 88,200
Variables: 16
$ year_id        <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2…
$ location_name  <chr> "Global", "Global", "Global", "Global", "Global",…
$ food           <chr> "fruits", "fruits", "fruits", "fruits", "fruits",…
$ age_group_id   <dbl> 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 30, 3…
$ age_group_name <chr> "25 to 29", "30 to 34", "35 to 39", "40 to 44", "…
$ sex            <chr> "Male", "Male", "Male", "Male", "Male", "Male", "…
$ parameter      <chr> "continuous", "continuous", "continuous", "contin…
$ mean           <dbl> 68.54567, 72.62209, 76.73157, 81.03910, 88.79202,…
$ upper          <dbl> 75.16806, 80.00433, 84.73205, 89.52747, 100.94301…
$ lower          <dbl> 63.17081, 66.78018, 71.04691, 73.82734, 80.40971,…
$ unit           <chr> "g/day", "g/day", "g/day", "g/day", "g/day", "g/d…
$ direction      <chr> "low", "low", "low", "low", "low", "low", "low", …
$ lower_optimal  <dbl> 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,…
$ optimal        <dbl> 250, 250, 250, 250, 250, 250, 250, 250, 250, 250,…
$ upper_optimal  <dbl> 300, 300, 300, 300, 300, 300, 300, 300, 300, 300,…
$ unit_optimal   <chr> "g", "g", "g", "g", "g", "g", "g", "g", "g", "g",…

It’s always a good idea to check that the values are what you expect after merging.

diet_and_guidelines %>%
  count(food, optimal)

# A tibble: 15 x 3
   food                        optimal     n
   <chr>                         <dbl> <int>
 1 calcium                        1.25   392
 2 fiber                         24      392
 3 fruits                       250      392
 4 legumes                       60      392
 5 milk                         435      392
 6 nuts and seeds                21      392
 7 polyunsaturated fatty acids   11      392
 8 processed meat                 2      392
 9 red meat                      23      392
10 seafood omega-3 fatty acids    0.25   392
11 sodium                         3      392
12 sugar-sweetened beverages      3      392
13 trans fatty acids              0.5    392
14 vegetables                   360      392
15 whole grains                 125      392

all_age_diet_and_guidelines %>%
  count(food, optimal)

# A tibble: 15 x 3
   food                        optimal     n
   <chr>                         <dbl> <int>
 1 calcium                        1.25  5880
 2 fiber                         24     5880
 3 fruits                       250     5880
 4 legumes                       60     5880
 5 milk                         435     5880
 6 nuts and seeds                21     5880
 7 polyunsaturated fatty acids   11     5880
 8 processed meat                 2     5880
 9 red meat                      23     5880
10 seafood omega-3 fatty acids    0.25  5880
11 sodium                         3     5880
12 sugar-sweetened beverages      3     5880
13 trans fatty acids              0.5   5880
14 vegetables                   360     5880
15 whole grains                 125     5880

guidelines

# A tibble: 15 x 6
   direction food          lower_optimal optimal upper_optimal unit_optimal
   <chr>     <chr>                 <dbl>   <dbl>         <dbl> <chr>       
 1 low       fruits                200    250            300   g           
 2 low       vegetables            290    360            430   g           
 3 low       legumes                50     60             70   g           
 4 low       whole grains          100    125            150   g           
 5 low       nuts and see…          16     21             25   g           
 6 low       milk                  350    435            520   g           
 7 high      red meat               18     23             27   g           
 8 high      processed me…           0      2              4   g           
 9 high      sugar-sweete…           0      3              5   g           
10 low       fiber                  19     24             28   g           
11 low       calcium                 1      1.25           1.5 g           
12 low       seafood omeg…           0.2    0.25           0.3 g           
13 low       polyunsatura…           9     11             13   %           
14 high      trans fatty …           0      0.5            1   %           
15 high      sodium                  1      3              5   g           

Looks good!

Calculate relative consumption

Again, we would like to compare the consumption rates of this dietary factors by different groups of people, but ideally we want to know this relative to the optimal guidelines.

Thus lets calcuate values of consuption that are relative to the suggested guidelines.

There are a few ways we could do this. One is to calculate a percentage of consumption based on the mean value for each observed value relative to the optimal value. To do this we will use the mutate() function of the dplyrpackage. This will create a new variable called mean_percent that will be equal to the division result of the mean variable value and the optimal variable multiplied by 100 to create a percentage.

diet_and_guidelines %<>%
  mutate(mean_percent = (mean/optimal)*100)

all_age_diet_and_guidelines %<>%
 mutate(mean_percent = (mean/optimal)*100)

Another option is to incorperate the range of optimal intakes and the direction that is associated with health risk. If the direction of risk is high and the consumption was greater than the optimal mean value, than the percentage is calculated based on the upper_optimal value, while if the direction of risk is low and the consumption is less than the optimal mean value, then the percentage is calculated based on the lower_optimal value. We will use the case_when() function of the dplyr package to do this. This allows us to specify values (indicated on the right side of the ~symbol) based on specific conditions (indicated on the left side of the ~ symbol). We can specify multiple conditions using the & symbol.

diet_and_guidelines %<>%
  mutate(range_percent =case_when(
  direction == "high" ~  (mean/upper_optimal)*100,
  direction == "low"  ~  (mean/lower_optimal)*100))

all_age_diet_and_guidelines %<>%
  mutate(range_percent =case_when(
  direction == "high" ~  (mean/upper_optimal)*100,
  direction == "low"  ~  (mean/lower_optimal)*100))


diet_and_guidelines %<>%
  mutate(percent_over_under =case_when(
  direction == "high" & mean>upper_optimal ~  ((mean-upper_optimal)/upper_optimal)*100,
  direction == "high" & mean<=upper_optimal ~ 0,
  direction == "low"  & mean>=lower_optimal ~ 0,
  direction == "low"  & mean<lower_optimal ~ ((lower_optimal-mean)/lower_optimal)*100))


all_age_diet_and_guidelines %<>%
  mutate(percent_over_under =case_when(
  direction == "high" & mean>upper_optimal ~  ((mean-upper_optimal)/upper_optimal)*100,
  direction == "high" & mean<=upper_optimal ~ 0,
  direction == "low"  & mean>=lower_optimal ~ 0,
  direction == "low"  & mean<lower_optimal ~ ((lower_optimal-mean)/lower_optimal)*100))

Yet another option is to create a binary outcome of if optimal consumption was achieved or not.

diet_and_guidelines %<>%
  mutate(opt_achieved = if_else(
    condition = direction =="low" & mean > lower_optimal |
                direction == "high" & mean < upper_optimal, 
    true = "Yes",
    false = "No"))

all_age_diet_and_guidelines %<>%
  mutate(opt_achieved = if_else(
    condition = direction =="low" & mean > lower_optimal |
                direction == "high" & mean < upper_optimal, 
    true = "Yes",
    false = "No"))

glimpse(diet_and_guidelines)

Observations: 5,880
Variables: 20
$ year_id            <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 201…
$ location_name      <chr> "Global", "Global", "Global", "Global", "Glob…
$ rei_id             <dbl> 111, 111, 112, 112, 113, 113, 114, 114, 115, …
$ food               <chr> "fruits", "fruits", "vegetables", "vegetables…
$ age_group_name     <chr> "All Available Ages", "All Available Ages", "…
$ sex                <chr> "Male", "Female", "Male", "Female", "Male", "…
$ parameter          <chr> "continuous", "continuous", "continuous", "co…
$ mean               <dbl> 8.703202e+01, 9.953182e+01, 1.975187e+02, 1.8…
$ upper              <dbl> 8.988734e+01, 1.023450e+02, 2.083679e+02, 1.9…
$ lower              <dbl> 8.452347e+01, 9.672140e+01, 1.887154e+02, 1.7…
$ unit               <chr> "g/day", "g/day", "g/day", "g/day", "g/day", …
$ direction          <chr> "low", "low", "low", "low", "low", "low", "lo…
$ lower_optimal      <dbl> 200.0, 200.0, 290.0, 290.0, 100.0, 100.0, 16.…
$ optimal            <dbl> 250.00, 250.00, 360.00, 360.00, 125.00, 125.0…
$ upper_optimal      <dbl> 300.0, 300.0, 430.0, 430.0, 150.0, 150.0, 25.…
$ unit_optimal       <chr> "g", "g", "g", "g", "g", "g", "g", "g", "g", …
$ mean_percent       <dbl> 34.8128083, 39.8127275, 54.8663179, 50.898279…
$ range_percent      <dbl> 43.5160103, 49.7659094, 68.1099118, 63.184071…
$ percent_over_under <dbl> 56.483990, 50.234091, 31.890088, 36.815929, 7…
$ opt_achieved       <chr> "No", "No", "No", "No", "No", "No", "No", "No…

glimpse(all_age_diet_and_guidelines)

Observations: 88,200
Variables: 20
$ year_id            <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 201…
$ location_name      <chr> "Global", "Global", "Global", "Global", "Glob…
$ food               <chr> "fruits", "fruits", "fruits", "fruits", "frui…
$ age_group_id       <dbl> 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 3…
$ age_group_name     <chr> "25 to 29", "30 to 34", "35 to 39", "40 to 44…
$ sex                <chr> "Male", "Male", "Male", "Male", "Male", "Male…
$ parameter          <chr> "continuous", "continuous", "continuous", "co…
$ mean               <dbl> 68.54567, 72.62209, 76.73157, 81.03910, 88.79…
$ upper              <dbl> 75.16806, 80.00433, 84.73205, 89.52747, 100.9…
$ lower              <dbl> 63.17081, 66.78018, 71.04691, 73.82734, 80.40…
$ unit               <chr> "g/day", "g/day", "g/day", "g/day", "g/day", …
$ direction          <chr> "low", "low", "low", "low", "low", "low", "lo…
$ lower_optimal      <dbl> 200, 200, 200, 200, 200, 200, 200, 200, 200, …
$ optimal            <dbl> 250, 250, 250, 250, 250, 250, 250, 250, 250, …
$ upper_optimal      <dbl> 300, 300, 300, 300, 300, 300, 300, 300, 300, …
$ unit_optimal       <chr> "g", "g", "g", "g", "g", "g", "g", "g", "g", …
$ mean_percent       <dbl> 27.41827, 29.04883, 30.69263, 32.41564, 35.51…
$ range_percent      <dbl> 34.27284, 36.31104, 38.36578, 40.51955, 44.39…
$ percent_over_under <dbl> 65.72716, 63.68896, 61.63422, 59.48045, 55.60…
$ opt_achieved       <chr> "No", "No", "No", "No", "No", "No", "No", "No…

One last thing that can be useful with data wrangling is to reshape the data into what is called the long format. This is very useful for creating visualizations with a very useful package called ggplot2.

To coerce an object into long format, we create more rows and fewer columns. For a more information about this, please see this case study.

We would like to put toether the different types of percentages of the optimal intake that we just calculated.

To get our data in long format we can use the pivot_longer() function of the dplyr package. We will also show how this would be done with the older version of this function, called gather().

For pivot_longer(), we will list the columns that we want to come together into the longer format using the cols argument. For gather() we would simply list the variables that we wish to consolidate. The names_to argument indicates the name of the variable that will include the character information about the values that we are consolidating, this is the variable names of the columns that we are bringing together. This is equivalent to the key aregument in gather(). The values_to is the name of the column that will contain the values of teh columns we are consolidating. This is equivalent to the value aregument in gather(). We can use contains() of the tidyr package to look at the variables with names that contain "percent" .

We would get an identical output from the methods.

diet_and_guidelines_long<-diet_and_guidelines %>%
pivot_longer(cols = contains("percent"), 
             names_to = "percent_type", 
             values_to = "percent")

diet_and_guidelines_long2<-diet_and_guidelines %>%
gather(contains("percent"),
       key = percent_type, 
       value = percent)

setequal(diet_and_guidelines_long, diet_and_guidelines_long2)

[1] TRUE

Let’s do the same for the age separated data.

all_age_diet_and_guidelines_long<-all_age_diet_and_guidelines %>%
pivot_longer(cols = contains("percent"), 
             names_to = "percent_type", 
             values_to = "percent")

Data Exploration

Exploring age collapsed data

Let’s take a look at the percent of consumption. Again we will use the base R summary() function:

diet_and_guidelines %>%
  select(mean_percent)%>%
  summary()

  mean_percent    
 Min.   :   0.02  
 1st Qu.:  11.20  
 Median :  34.71  
 Mean   : 187.86  
 3rd Qu.:  73.12  
 Max.   :7693.68  

Wow! Some of the values are nearly zero, suggesting that some people are consuming basically zero percent of what is suggested for optimal health. On the other hand, for some dietary factors people are consuming over 7,000 percent what is suggested!

This is why it is important to look at the direction of consumption that could be harmful. For example if there is a population that consumes large amounts of vegatables this could be a good thing, but if there is a population consuming large amounts of sodium this would be a bad thing.

Let’s take a look to see what dietary factors are at the extremes by arranging the data using the arrange() function of the dplyr package. We can arrange by smallest to largest using the default and we can arrange largest to smallest using the minus sign -.

diet_and_guidelines %>%
  arrange(-mean_percent)%>%
  glimpse()

Observations: 5,880
Variables: 20
$ year_id            <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 201…
$ location_name      <chr> "Costa Rica", "Colombia", "Panama", "Venezuel…
$ rei_id             <dbl> 118, 118, 118, 118, 118, 118, 118, 118, 118, …
$ food               <chr> "sugar-sweetened beverages", "sugar-sweetened…
$ age_group_name     <chr> "All Available Ages", "All Available Ages", "…
$ sex                <chr> "Male", "Male", "Male", "Male", "Male", "Male…
$ parameter          <chr> "continuous", "continuous", "continuous", "co…
$ mean               <dbl> 230.8105, 201.7180, 194.4116, 193.9048, 188.4…
$ upper              <dbl> 247.9342, 216.9260, 208.8040, 208.8702, 203.3…
$ lower              <dbl> 215.5919, 189.4812, 181.7150, 180.7209, 174.5…
$ unit               <chr> "g/day", "g/day", "g/day", "g/day", "g/day", …
$ direction          <chr> "high", "high", "high", "high", "high", "high…
$ lower_optimal      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ optimal            <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, …
$ upper_optimal      <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
$ unit_optimal       <chr> "g", "g", "g", "g", "g", "g", "g", "g", "g", …
$ mean_percent       <dbl> 7693.684, 6723.934, 6480.386, 6463.492, 6282.…
$ range_percent      <dbl> 4616.211, 4034.361, 3888.232, 3878.095, 3769.…
$ percent_over_under <dbl> 4516.211, 3934.361, 3788.232, 3778.095, 3669.…
$ opt_achieved       <chr> "No", "No", "No", "No", "No", "No", "No", "No…

Ok, so it looks like sugar-sweetened beverages are really overconsumed in some parts of the world!

Recall from the supplementary table from the article that overconsumption of sugar-sweetened beverages is associated with both Diabetes mellitus type 2 and Ischemic heart disease. This article discusses some of the contrevsy over the potential health risks associated with high consumption of sugar.

It still looks quite bad if we look at the other calculated percentage values.

diet_and_guidelines %>%
   select(contains("percent"))%>%
   summary()

  mean_percent     range_percent      percent_over_under
 Min.   :   0.02   Min.   :   0.025   Min.   :   0.00   
 1st Qu.:  11.20   1st Qu.:  13.263   1st Qu.:   0.00   
 Median :  34.71   Median :  38.941   Median :  56.82   
 Mean   : 187.86   Mean   : 127.272   Mean   : 122.91   
 3rd Qu.:  73.12   3rd Qu.:  72.414   3rd Qu.:  87.78   
 Max.   :7693.68   Max.   :4616.211   Max.   :4516.21   

So some places are still consuming 4,000 percent more than the upper range of the suggested optimal intake.

Let’s take a look at global levels:

diet_and_guidelines %>%
  filter(food =="sugar-sweetened beverages" & location_name =="Global")

# A tibble: 2 x 20
  year_id location_name rei_id food  age_group_name sex   parameter  mean
    <dbl> <chr>          <dbl> <chr> <chr>          <chr> <chr>     <dbl>
1    2017 Global           118 suga… All Available… Male  continuo…  65.5
2    2017 Global           118 suga… All Available… Fema… continuo…  47.7
# … with 12 more variables: upper <dbl>, lower <dbl>, unit <chr>,
#   direction <chr>, lower_optimal <dbl>, optimal <dbl>,
#   upper_optimal <dbl>, unit_optimal <chr>, mean_percent <dbl>,
#   range_percent <dbl>, percent_over_under <dbl>, opt_achieved <chr>

For those who are less familiar with the metric system where grams are equivalent to milliliters, it may be useful to realize how many fluid ounces the max amount of consumption per day (~248) grams actaully is.

There are 0.35247 ounces in one gram.

#top amount in ounces
0.35247*247.9342 

[1] 87.38937

Ok, so the top consumers are drinking about 87 fluid ounces per day. Since there are 12 ounces in a single can of soda, this is about 7.25 sodas per day. Globally on average, males are drinking 1.9238988 about sodas worth of sweetened beverages, while females are drinking about 1.4010682.

Let’s take a look at what is underconsumed:

diet_and_guidelines %>%
  arrange(mean_percent)%>%
  glimpse()

Observations: 5,880
Variables: 20
$ year_id            <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 201…
$ location_name      <chr> "Chad", "Chad", "Mali", "Mali", "Burkina Faso…
$ rei_id             <dbl> 122, 122, 122, 122, 122, 122, 122, 122, 122, …
$ food               <chr> "polyunsaturated fatty acids", "polyunsaturat…
$ age_group_name     <chr> "All Available Ages", "All Available Ages", "…
$ sex                <chr> "Female", "Male", "Female", "Male", "Female",…
$ parameter          <chr> "continuous", "continuous", "continuous", "co…
$ mean               <dbl> 0.002227074, 0.002295916, 0.002301266, 0.0023…
$ upper              <dbl> 0.002383410, 0.002439931, 0.002461209, 0.0025…
$ lower              <dbl> 0.002082132, 0.002161922, 0.002152212, 0.0022…
$ unit               <chr> "%energy/day", "%energy/day", "%energy/day", …
$ direction          <chr> "low", "low", "low", "low", "low", "low", "lo…
$ lower_optimal      <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, …
$ optimal            <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 1…
$ upper_optimal      <dbl> 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 1…
$ unit_optimal       <chr> "%", "%", "%", "%", "%", "%", "%", "%", "%", …
$ mean_percent       <dbl> 0.02024613, 0.02087196, 0.02092060, 0.0215757…
$ range_percent      <dbl> 0.02474527, 0.02551018, 0.02556962, 0.0263703…
$ percent_over_under <dbl> 99.97525, 99.97449, 99.97443, 99.97363, 99.96…
$ opt_achieved       <chr> "No", "No", "No", "No", "No", "No", "No", "No…

#glimpse(filter(diet_and_guidelines, location_name == "Laos", mean>200))

On the otherhand, it looks like some places are consuming almost no polyunsaturated fatty acids. These are fats that found in plant-based sources like seeds and nuts. According to an article about polyunsaturated fatty acids and its influence on health:

Coronary heart disease (CHD) is the leading cause of death worldwide … The types of dietary fats consumed play an important role in CHD risk, representing key modifiable risk factors…In particular, higher intakes of trans fat (TFA) and of saturated fat (SFA) replacing ω‐6 (n‐6) polyunsaturated fat (PUFA) are associated with increased CHD… whereas higher intake of PUFA replacing either SFA or carbohydrate is associated with lower risk.

Let’s get an idea about how countries compare in terms of how many of the dietary factors are consumed at the optimal level (the opt_achieved variable).

diet_and_guidelines %>%
  count(opt_achieved)

# A tibble: 2 x 2
  opt_achieved     n
  <chr>        <int>
1 No            4360
2 Yes           1520

Looks like overal, only 34.8623853 of dietary factors for all tested populations were at optimal levels.

Let’s get an idea about how countries compare.

diet_and_guidelines %>%
  count(opt_achieved, location_name) %>%
  filter(opt_achieved == "Yes") %>%
  arrange(-n) %>%
  print(n = 1e3)

# A tibble: 196 x 3
    opt_achieved location_name                        n
    <chr>        <chr>                            <int>
  1 Yes          Qatar                               13
  2 Yes          Rwanda                              13
  3 Yes          Turkey                              13
  4 Yes          Dominican Republic                  12
  5 Yes          Egypt                               12
  6 Yes          Malawi                              12
  7 Yes          Benin                               11
  8 Yes          Burkina Faso                        11
  9 Yes          Cameroon                            11
 10 Yes          Cuba                                11
 11 Yes          Ethiopia                            11
 12 Yes          Ghana                               11
 13 Yes          Lesotho                             11
 14 Yes          Mali                                11
 15 Yes          Nigeria                             11
 16 Yes          Senegal                             11
 17 Yes          Togo                                11
 18 Yes          Tunisia                             11
 19 Yes          Angola                              10
 20 Yes          Azerbaijan                          10
 21 Yes          Bahrain                             10
 22 Yes          Barbados                            10
 23 Yes          Bermuda                             10
 24 Yes          Chad                                10
 25 Yes          Comoros                             10
 26 Yes          Cote d'Ivoire                       10
 27 Yes          El Salvador                         10
 28 Yes          Guatemala                           10
 29 Yes          Haiti                               10
 30 Yes          India                               10
 31 Yes          Iran                                10
 32 Yes          Kenya                               10
 33 Yes          Lebanon                             10
 34 Yes          Mauritania                          10
 35 Yes          Mauritius                           10
 36 Yes          Morocco                             10
 37 Yes          Mozambique                          10
 38 Yes          Myanmar                             10
 39 Yes          Nepal                               10
 40 Yes          Nicaragua                           10
 41 Yes          Niger                               10
 42 Yes          North Korea                         10
 43 Yes          Sao Tome and Principe               10
 44 Yes          Sierra Leone                        10
 45 Yes          South Sudan                         10
 46 Yes          Tanzania                            10
 47 Yes          Timor-Leste                         10
 48 Yes          Uganda                              10
 49 Yes          United Arab Emirates                10
 50 Yes          Zambia                              10
 51 Yes          Armenia                              9
 52 Yes          Bosnia and Herzegovina               9
 53 Yes          Burundi                              9
 54 Yes          Costa Rica                           9
 55 Yes          Equatorial Guinea                    9
 56 Yes          Eritrea                              9
 57 Yes          Fiji                                 9
 58 Yes          Honduras                             9
 59 Yes          Iraq                                 9
 60 Yes          Kuwait                               9
 61 Yes          Oman                                 9
 62 Yes          Papua New Guinea                     9
 63 Yes          Solomon Islands                      9
 64 Yes          Tajikistan                           9
 65 Yes          Trinidad and Tobago                  9
 66 Yes          Uzbekistan                           9
 67 Yes          Virgin Islands, U.S.                 9
 68 Yes          Afghanistan                          8
 69 Yes          Albania                              8
 70 Yes          Algeria                              8
 71 Yes          American Samoa                       8
 72 Yes          Antigua and Barbuda                  8
 73 Yes          Bangladesh                           8
 74 Yes          Belize                               8
 75 Yes          Bhutan                               8
 76 Yes          Botswana                             8
 77 Yes          Brazil                               8
 78 Yes          Cambodia                             8
 79 Yes          Cape Verde                           8
 80 Yes          Congo                                8
 81 Yes          Democratic Republic of the Congo     8
 82 Yes          Djibouti                             8
 83 Yes          Dominica                             8
 84 Yes          Federated States of Micronesia       8
 85 Yes          Gabon                                8
 86 Yes          Georgia                              8
 87 Yes          Grenada                              8
 88 Yes          Guam                                 8
 89 Yes          Guinea                               8
 90 Yes          Guinea-Bissau                        8
 91 Yes          Guyana                               8
 92 Yes          Indonesia                            8
 93 Yes          Israel                               8
 94 Yes          Jamaica                              8
 95 Yes          Jordan                               8
 96 Yes          Kazakhstan                           8
 97 Yes          Kiribati                             8
 98 Yes          Kyrgyzstan                           8
 99 Yes          Laos                                 8
100 Yes          Liberia                              8
101 Yes          Libya                                8
102 Yes          Macedonia                            8
103 Yes          Madagascar                           8
104 Yes          Malaysia                             8
105 Yes          Maldives                             8
106 Yes          Marshall Islands                     8
107 Yes          Mexico                               8
108 Yes          Moldova                              8
109 Yes          Namibia                              8
110 Yes          Northern Mariana Islands             8
111 Yes          Pakistan                             8
112 Yes          Palestine                            8
113 Yes          Peru                                 8
114 Yes          Philippines                          8
115 Yes          Saint Lucia                          8
116 Yes          Saint Vincent and the Grenadines     8
117 Yes          Samoa                                8
118 Yes          Saudi Arabia                         8
119 Yes          Seychelles                           8
120 Yes          Somalia                              8
121 Yes          Sri Lanka                            8
122 Yes          Sudan                                8
123 Yes          Suriname                             8
124 Yes          Swaziland                            8
125 Yes          Syria                                8
126 Yes          Taiwan                               8
127 Yes          Thailand                             8
128 Yes          The Bahamas                          8
129 Yes          The Gambia                           8
130 Yes          Tonga                                8
131 Yes          Turkmenistan                         8
132 Yes          Vanuatu                              8
133 Yes          Yemen                                8
134 Yes          Zimbabwe                             8
135 Yes          Bolivia                              7
136 Yes          Central African Republic             7
137 Yes          Luxembourg                           7
138 Yes          Montenegro                           7
139 Yes          Panama                               7
140 Yes          Paraguay                             7
141 Yes          Puerto Rico                          7
142 Yes          South Africa                         7
143 Yes          Ukraine                              7
144 Yes          Venezuela                            7
145 Yes          Vietnam                              7
146 Yes          Canada                               6
147 Yes          Colombia                             6
148 Yes          Greece                               6
149 Yes          Japan                                6
150 Yes          Malta                                6
151 Yes          Mongolia                             6
152 Yes          Portugal                             6
153 Yes          Romania                              6
154 Yes          Serbia                               6
155 Yes          Argentina                            5
156 Yes          Belarus                              5
157 Yes          Bulgaria                             5
158 Yes          Croatia                              5
159 Yes          Ecuador                              5
160 Yes          Finland                              5
161 Yes          Global                               5
162 Yes          Iceland                              5
163 Yes          Netherlands                          5
164 Yes          Russian Federation                   5
165 Yes          Sweden                               5
166 Yes          Andorra                              4
167 Yes          Australia                            4
168 Yes          Austria                              4
169 Yes          Belgium                              4
170 Yes          Brunei                               4
171 Yes          Chile                                4
172 Yes          China                                4
173 Yes          Cyprus                               4
174 Yes          Denmark                              4
175 Yes          Estonia                              4
176 Yes          France                               4
177 Yes          Germany                              4
178 Yes          Ireland                              4
179 Yes          Italy                                4
180 Yes          Latvia                               4
181 Yes          Lithuania                            4
182 Yes          New Zealand                          4
183 Yes          Norway                               4
184 Yes          Poland                               4
185 Yes          Singapore                            4
186 Yes          South Korea                          4
187 Yes          Spain                                4
188 Yes          Switzerland                          4
189 Yes          United Kingdom                       4
190 Yes          Uruguay                              4
191 Yes          Czech Republic                       3
192 Yes          Greenland                            3
193 Yes          Hungary                              3
194 Yes          Slovakia                             3
195 Yes          Slovenia                             3
196 Yes          United States                        3

Looks as though on average the populations (both male and female separately) in Qatar, Rwanda, and Turkey consumped the optimal level of intake for the largest number of dietary factors (13 out of 30 (for the 15 dietary factors for males and females)).

In contrast, the Czech Republic, Greenland, Hungary, Slovakia, Slovenia, and the United States had the poorest consumption rates (27 out of 30 were not at optimal levels).

diet_and_guidelines %>%
  count(opt_achieved, location_name) %>%
  filter(opt_achieved == "No") %>%
  arrange(-n) %>%
  print(n = 1e3)

# A tibble: 196 x 3
    opt_achieved location_name                        n
    <chr>        <chr>                            <int>
  1 No           Czech Republic                      27
  2 No           Greenland                           27
  3 No           Hungary                             27
  4 No           Slovakia                            27
  5 No           Slovenia                            27
  6 No           United States                       27
  7 No           Andorra                             26
  8 No           Australia                           26
  9 No           Austria                             26
 10 No           Belgium                             26
 11 No           Brunei                              26
 12 No           Chile                               26
 13 No           China                               26
 14 No           Cyprus                              26
 15 No           Denmark                             26
 16 No           Estonia                             26
 17 No           France                              26
 18 No           Germany                             26
 19 No           Ireland                             26
 20 No           Italy                               26
 21 No           Latvia                              26
 22 No           Lithuania                           26
 23 No           New Zealand                         26
 24 No           Norway                              26
 25 No           Poland                              26
 26 No           Singapore                           26
 27 No           South Korea                         26
 28 No           Spain                               26
 29 No           Switzerland                         26
 30 No           United Kingdom                      26
 31 No           Uruguay                             26
 32 No           Argentina                           25
 33 No           Belarus                             25
 34 No           Bulgaria                            25
 35 No           Croatia                             25
 36 No           Ecuador                             25
 37 No           Finland                             25
 38 No           Global                              25
 39 No           Iceland                             25
 40 No           Netherlands                         25
 41 No           Russian Federation                  25
 42 No           Sweden                              25
 43 No           Canada                              24
 44 No           Colombia                            24
 45 No           Greece                              24
 46 No           Japan                               24
 47 No           Malta                               24
 48 No           Mongolia                            24
 49 No           Portugal                            24
 50 No           Romania                             24
 51 No           Serbia                              24
 52 No           Bolivia                             23
 53 No           Central African Republic            23
 54 No           Luxembourg                          23
 55 No           Montenegro                          23
 56 No           Panama                              23
 57 No           Paraguay                            23
 58 No           Puerto Rico                         23
 59 No           South Africa                        23
 60 No           Ukraine                             23
 61 No           Venezuela                           23
 62 No           Vietnam                             23
 63 No           Afghanistan                         22
 64 No           Albania                             22
 65 No           Algeria                             22
 66 No           American Samoa                      22
 67 No           Antigua and Barbuda                 22
 68 No           Bangladesh                          22
 69 No           Belize                              22
 70 No           Bhutan                              22
 71 No           Botswana                            22
 72 No           Brazil                              22
 73 No           Cambodia                            22
 74 No           Cape Verde                          22
 75 No           Congo                               22
 76 No           Democratic Republic of the Congo    22
 77 No           Djibouti                            22
 78 No           Dominica                            22
 79 No           Federated States of Micronesia      22
 80 No           Gabon                               22
 81 No           Georgia                             22
 82 No           Grenada                             22
 83 No           Guam                                22
 84 No           Guinea                              22
 85 No           Guinea-Bissau                       22
 86 No           Guyana                              22
 87 No           Indonesia                           22
 88 No           Israel                              22
 89 No           Jamaica                             22
 90 No           Jordan                              22
 91 No           Kazakhstan                          22
 92 No           Kiribati                            22
 93 No           Kyrgyzstan                          22
 94 No           Laos                                22
 95 No           Liberia                             22
 96 No           Libya                               22
 97 No           Macedonia                           22
 98 No           Madagascar                          22
 99 No           Malaysia                            22
100 No           Maldives                            22
101 No           Marshall Islands                    22
102 No           Mexico                              22
103 No           Moldova                             22
104 No           Namibia                             22
105 No           Northern Mariana Islands            22
106 No           Pakistan                            22
107 No           Palestine                           22
108 No           Peru                                22
109 No           Philippines                         22
110 No           Saint Lucia                         22
111 No           Saint Vincent and the Grenadines    22
112 No           Samoa                               22
113 No           Saudi Arabia                        22
114 No           Seychelles                          22
115 No           Somalia                             22
116 No           Sri Lanka                           22
117 No           Sudan                               22
118 No           Suriname                            22
119 No           Swaziland                           22
120 No           Syria                               22
121 No           Taiwan                              22
122 No           Thailand                            22
123 No           The Bahamas                         22
124 No           The Gambia                          22
125 No           Tonga                               22
126 No           Turkmenistan                        22
127 No           Vanuatu                             22
128 No           Yemen                               22
129 No           Zimbabwe                            22
130 No           Armenia                             21
131 No           Bosnia and Herzegovina              21
132 No           Burundi                             21
133 No           Costa Rica                          21
134 No           Equatorial Guinea                   21
135 No           Eritrea                             21
136 No           Fiji                                21
137 No           Honduras                            21
138 No           Iraq                                21
139 No           Kuwait                              21
140 No           Oman                                21
141 No           Papua New Guinea                    21
142 No           Solomon Islands                     21
143 No           Tajikistan                          21
144 No           Trinidad and Tobago                 21
145 No           Uzbekistan                          21
146 No           Virgin Islands, U.S.                21
147 No           Angola                              20
148 No           Azerbaijan                          20
149 No           Bahrain                             20
150 No           Barbados                            20
151 No           Bermuda                             20
152 No           Chad                                20
153 No           Comoros                             20
154 No           Cote d'Ivoire                       20
155 No           El Salvador                         20
156 No           Guatemala                           20
157 No           Haiti                               20
158 No           India                               20
159 No           Iran                                20
160 No           Kenya                               20
161 No           Lebanon                             20
162 No           Mauritania                          20
163 No           Mauritius                           20
164 No           Morocco                             20
165 No           Mozambique                          20
166 No           Myanmar                             20
167 No           Nepal                               20
168 No           Nicaragua                           20
169 No           Niger                               20
170 No           North Korea                         20
171 No           Sao Tome and Principe               20
172 No           Sierra Leone                        20
173 No           South Sudan                         20
174 No           Tanzania                            20
175 No           Timor-Leste                         20
176 No           Uganda                              20
177 No           United Arab Emirates                20
178 No           Zambia                              20
179 No           Benin                               19
180 No           Burkina Faso                        19
181 No           Cameroon                            19
182 No           Cuba                                19
183 No           Ethiopia                            19
184 No           Ghana                               19
185 No           Lesotho                             19
186 No           Mali                                19
187 No           Nigeria                             19
188 No           Senegal                             19
189 No           Togo                                19
190 No           Tunisia                             19
191 No           Dominican Republic                  18
192 No           Egypt                               18
193 No           Malawi                              18
194 No           Qatar                               17
195 No           Rwanda                              17
196 No           Turkey                              17

Let’s look at the raw US data:

diet_and_guidelines %>%
  filter(location_name == "United States") %>%
  glimpse()

Observations: 30
Variables: 20
$ year_id            <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 2017, 201…
$ location_name      <chr> "United States", "United States", "United Sta…
$ rei_id             <dbl> 111, 111, 112, 112, 113, 113, 114, 114, 115, …
$ food               <chr> "fruits", "fruits", "vegetables", "vegetables…
$ age_group_name     <chr> "All Available Ages", "All Available Ages", "…
$ sex                <chr> "Male", "Female", "Male", "Female", "Male", "…
$ parameter          <chr> "continuous", "continuous", "continuous", "co…
$ mean               <dbl> 114.24773940, 129.11513300, 208.93230110, 192…
$ upper              <dbl> 121.37292480, 136.94611510, 217.80223830, 200…
$ lower              <dbl> 107.89222320, 122.33941970, 200.05582000, 184…
$ unit               <chr> "g/day", "g/day", "g/day", "g/day", "g/day", …
$ direction          <chr> "low", "low", "low", "low", "low", "low", "lo…
$ lower_optimal      <dbl> 200.0, 200.0, 290.0, 290.0, 100.0, 100.0, 16.…
$ optimal            <dbl> 250.00, 250.00, 360.00, 360.00, 125.00, 125.0…
$ upper_optimal      <dbl> 300.0, 300.0, 430.0, 430.0, 150.0, 150.0, 25.…
$ unit_optimal       <chr> "g", "g", "g", "g", "g", "g", "g", "g", "g", …
$ mean_percent       <dbl> 45.699096, 51.646053, 58.036750, 53.515649, 2…
$ range_percent      <dbl> 57.123870, 64.557566, 72.045621, 66.433220, 2…
$ percent_over_under <dbl> 42.876130, 35.442434, 27.954379, 33.566780, 7…
$ opt_achieved       <chr> "No", "No", "No", "No", "No", "No", "No", "No…

Let’s see how males and females compare for achieving the optimal intake:

count(diet_and_guidelines, sex, opt_achieved)

# A tibble: 4 x 3
  sex    opt_achieved     n
  <chr>  <chr>        <int>
1 Female No            2171
2 Female Yes            769
3 Male   No            2189
4 Male   Yes            751

Looks pretty similar, but it may be a bit better for females. We will evaluate this further.

Here is a way we can visualise this with the ggplot2 package. The ggplot2 package creates plots by using layers. Notice in the following code how there is a plus sign between the ggplot() function and the geom_histogram() function. With ggplot2 we select what data we would like to plot using the first function (ggplot()) and then we add on additional layers of complexity (these layers can even involve different data). The aes() argument specifies what aspects of the data will be plotted where. the geom_* function specifies what type of plot to create (e.g. geom_histogram() create a histogram).

We will see later how we can add many layers to plots with ggplot2. For additional information on using ggplot2, see this case study.

diet_and_guidelines %>%
  ggplot(aes(opt_achieved , col= sex))+
  geom_bar()

Continuing with ggplot2 we will now create a different plot - this time we will create a series of boxplots. We will use the facet_wrap() function of ggplot2 to allow us to create many different plots simultaneously. In this case we can look at boxplots for the different dietary factors colored by sex. The scales argument when set to "free" means that each of the sequentual plot created by the facet can have a differnt scale for the y axis, otherwise, by default they are constrained to the same scale.

diet_and_guidelines%<>%
  mutate(food_to_plot =
  str_replace( 
  string =pull(diet_and_guidelines,food), 
  pattern = " ", 
  replacement = "\n"))

diet_and_guidelines %>%
  ggplot(aes(y = mean_percent , x= sex, color = sex))+
  geom_boxplot()+
  facet_wrap(~food_to_plot, scales = "free", nrow = 3, strip.position = "right")+
  theme(strip.text.y  =  element_text(size = 8),
        axis.text.x = element_text(angle = 70, hjust = 1))

If we just look at differences by sex for the specific dietary factors, males appear to potentially consume more of many of the factors, including possibly more sodium, fiber, calcium, red meat, and sugar-sweetened beverages than females. Females may consume more fruit.

Exploring the data separated by age

Now we will take a look at the data that is separated by age groups.

First, recall that we have 15 different age groups starting from age 25 to 95 plus.

all_age_diet_and_guidelines %>%
 count(age_group_name)

# A tibble: 15 x 2
   age_group_name     n
   <chr>          <int>
 1 25 to 29        5880
 2 30 to 34        5880
 3 35 to 39        5880
 4 40 to 44        5880
 5 45 to 49        5880
 6 50 to 54        5880
 7 55 to 59        5880
 8 60 to 64        5880
 9 65 to 69        5880
10 70 to 74        5880
11 75 to 79        5880
12 80 to 84        5880
13 85 to 89        5880
14 90 to 94        5880
15 95 plus         5880

sep_age_diet_data %>%
  ggplot(aes(y = mean , x= age_group_name, col = sex))+
  geom_boxplot()+
  facet_wrap(~food, scales = "free", nrow = 6)+
  theme(axis.text.x = element_text(angle = 70, hjust = 1))

# sep_age_diet_data %>%
#   ggplot(aes(y = mean , x= age_group_name, col = sex))+
#   geom_boxplot()+
#   facet_wrap(~food, scales = "free")+
#   facet_wrap(~food, scales = "free", nrow = 3, strip.position = "right")+
#   theme(strip.text.y  =  element_text(size = 8),
#         axis.text.x = element_text(angle = 70, hjust = 1))

We can see from these plots that there appears to be age differences and gender differences for some of the different dietary factors. We will work to create clearer figures later on. However these figures have given us a better sense of the data that we are working with.

Data Analysis

Recall what our main question were:

Our main questions are:

1. What are the global trends for potentially harmful diets?
2. How do males and females compare?
3. How do different age groups compare for these dietary factors?
4. How do different countries compare? In particular, how does the US compare to other contries in terms of diet trends?

We have some general sense about global trends for the risk-associated dietary factors, however we want to know more.

We are interested in how the 2 genders compare, how the 195 different countries compare and how the 15 differnt age groups compare. We may have learned that we can compare two groups using a \(t\)-test, but how can we compare more than 2 groups? For more information on the \(t\)-test see this case study.

In order to make inference about these comparisons, it is helpful to perform statistical tests. First we are going to talk about a statistical method called regression.

Perhaps you have heard about the ANOVA test which can be useful for testing more than two groups. ANOVA stands for “ANalysis Of VAriance”. It, just like the \(t\)-test, is a specialized type of regression. More on that to come…

Regression

So what is regression?

The term was coined in 1877 by Sir Francis Galton (Charles Darwin’s half-cousin!) in his article about the relationship between heriditary traits and population averages. He particularly focused on height and kinship or relatedness. The word in general means to go back to a simpler mode. Galton noticed that individuals with parents who had an extreme trait, such as height, tended to have a height more similar to the average of the population rather than the extreme height of their parents. For example if parents were very tall, their children were likely to be a bit shorter than their parents and therefore closer to the population average. Thus the children regressed towards the mean or in Galton’s words the offspring showed:

“a regression towards mediocrity”

See here for more information about this history.

When we think about this from a statistical standpoint, regression allows us to estimate or regress relationships between variables to a “simple” model. We do this by estimating the mean of an outcome. This can be useful for predicting future values of the outcome based on the approximation of the real relationship between the varaibles within the model.

Using the ordinary least squares method we can identify a line that best fits the data by minimizing the sum of the squared distances between each point and the line.

Fitting a line to the data like this allows us to create a formula for the line using an intercept and a slope, so that we can then estimate mean values of \(Y\) (dependent/outcome variable) given known values of \(X\) (independent/predictor/covariate/explanatory variable(s)). People will also say that we are "regressing \(Y\) on \(X\).

You may have seen the formula for a line written like this:

\[Y = mX + b\]

or

\[Y = aX + b\]

In this case \(m\) or \(a\) is the slope of the line and \(b\) is a constant and represents the intercept or the point where the y axis is crossed by the line, when \(x = 0\).

We can also write this model like this:

\[Y = \beta_{1}X +\beta_{0}\]

Now \(\beta_{1}\) called Beta one is our slope notation and \(B0\) Beta not is our intercept notation.

Importantly the slope (\(m\) or \(a\) or \(\beta_{1}\)) gives us a measure of the strength of the influence of the independent variable(s)(\(X_{i}\)) on the dependent variable (\(Y\)).

Check out this interactive explanation of how the ordinary least squares method works.

Here is an image of what we are saying about the ordinary least squares regression to fit a line to data:

[source]

In some cases we can fit a line perfectly and all points will lie on the line with no distance to the line:

In this case, the slope or \(\beta_{1}\) is 1 and the intercept is \(\beta_{0}\) 10. We can see that if \(X\) were 50, \(Y\) would be approximately 60. This is very unusal in statistical analysis however, as often the relationship between variables is more complicated.

In other cases there will be greater distances between the line and the points. Like this regression:

This concept has been extended to allow for comparisons of different types of outcomes (the ys), to involve various numbers of covariates (indpendent variables - the xs), and to allow for different shapes of lines. For a guide on how to perform regressions in R see here.

So let’s get back to our data…What types of outcomes do we have?

Our outcomes (\(Y\)):

We can evaluate the raw consumption or the percent of optimal consumption values that we calculated. These outcomes would all be what we call continuous beause our values can take on any numeric value within the range of possible values. Binary outcomes in contrast, have only two possible categorical values. We can also evaluate our binary outcome of if the optimal level of consumption was achieved or not.

Continuous outcomes can be evaluated with linear regression, while binary outcomes can be evaluated with logistic regression. See here for a guide on different types of regression methods.

We have both continuous outcomes of raw consumption rates or percentages of consumption relatative to the optimal guideline amounts as well as binary outcomes of “yes” or “no” about if a population achieved the optimal amount of consumption.

Our covariates (\(X_{i}\)):

We have data about age group, sex, and region informationfor all of our outcomes.

Again we are intereseted in comparing the sexes, the different age groups, and the different countries for their consumption of the different dietary factors. We can apply regression to determine if there is any influence of these group identities on the dietary consumption outcomes.

Before we get started, let’s remove the global values from our data and set them aside, as this is really a mean of all the country values.

global <- diet_and_guidelines %>% filter(location_name == "Global")
diet_and_guidelines %<>% filter(location_name != "Global")

Influence of sex on dietary outcomes

Let’s focus on the dietary factors that appeared to potentially have a difference between genders based on our figure in our exploratory analysis.

“If we just look at differences by sex for the specific dietary factors, males appear to potentially consume more of many of the factors, including possibly more sodium, fiber, calcium, red meat, and sugar-sweetened beverages than females. Females may consume more fruit.”

First let’s take a look at red meat.

We can compare the red meat consumption of males and females around the world using the well known \(t\)-test using the t.test() function and a linear regression model using the lm() function (both are included in stats package that is installed with R) and we will get the same results. See here for additional explanation about why that is the case. Here and here are also great sources about how many commonly known statistical tests are specialized forms of regression.

Before we get started, let’s think about the assumptions of both tests.

\(t\)-test assumptions:

1. Normality of the data for both groups (this is not as much of an issue if the number of observations is relatively large total n>30 - which is indeed the case for us!)
2. Equal variance between the two groups (make sure you do the correct test if the data is not normal)
3. Balanced sample sizes of the two groups (got that!)
4. Independent observations (or independent paired observations - got that too!)

We can evaluate if our data is normally distributed by plotting the distribution and by creating Q-Q plots.

If our data is not normally distributed, we can consider these options:

1. We can still perform a t-test if our n is large
2. We can transform the data before performing a t-test
3. We can use a nonparametric test (Wilcoxon signed rank test, the Wilcoxon rank sum test, and the Two-sample Kolmogorov-Smirnov (KS) test)
4. We can perform a t-test with resampling methods (which should be especially considered when the groups are imbalanced)

See this case study for more information on \(t\)-test assumptions.

Linear regression assumptions

L (linear) - There is a linear relationship between the variables. I (independent) - The samples are independent from one another. N (normal) - The residuals are normally distributed. E (equal variances) - The variance of the groups is similar.

Distributions

In order to apply a statistical test to compare the means, one of the first things to do is to explore the frequency of the different observed values. One way to summarize the frequency of different observed values is the frequency distribution, which can be shown in a table or a plot. See here for more information about distributions.

We will use the geom_histogram() of the ggplot2package to create a histogram to evaluate the frequency distributions of our data. The facet_wrap() function of the ggplot2 package allows us to look at different parts of our data in separate plots.

diet_and_guidelines %>%
  filter(food =="red meat") %>%
  ggplot(aes(x=mean_percent)) +
  geom_histogram() +
  facet_wrap(~ sex) 

Now we add another dimension to our facet with the plus sign. Now it becomes very useful that we reformmated our data to the long format, as it makes it easy to facet across the different percentages that we calculated. We can specify that we want to have different y axis scales using the scales argument as: scales = "free".

diet_and_guidelines_long %>%
  filter(food =="red meat") %>%
  ggplot(aes(x = percent)) +
  geom_histogram() +
  facet_wrap(~ percent_type + sex, scales = "free")

diet_and_guidelines_long %>%
  filter(food =="red meat") %>%
  ggplot(aes(sample = percent)) +
  facet_wrap(~ percent_type + sex)+
  geom_qq() +
  geom_qq_line()

Looks like our data is right skewed.

We can transform our data to make it more normally distrubuted. When data is highly right skewed, a log transformation is often helpful.

Let’s take a look a the log of our percent of optimal consumption values.

diet_and_guidelines_long %>%
    filter(food =="red meat") %>%
  ggplot(aes(x = log10(percent))) +
  geom_histogram() +
  facet_wrap(~ percent_type)

diet_and_guidelines_long %>%
    filter(food =="red meat") %>%
  ggplot(aes(sample = log10(percent))) +
  facet_wrap(~ percent_type)+
  geom_qq() +
  geom_qq_line()

How about the raw consumption data:

diet_and_guidelines_long %>%
  filter(food == "red meat") %>%
  ggplot(aes(x = (mean))) +
  geom_histogram()

Also right skewed… so let’s try taking the log of the mean value:

diet_and_guidelines_long %>%
  filter(food == "red meat") %>%
  ggplot(aes(x = (log10(mean)))) +
  geom_histogram()

This is better too.

Ok, so we can get our data to look fairly normal, which is good for our \(t\)-test assumptions. The other thing we need to check is if the variance in red meat consumption is similar between the two genders. We can use the var.test() of the stats package using the log normalized data, as this data is fairly normally distributed. We can use the mood.test on the raw data which is skewed.

mood.test(pull(filter(diet_and_guidelines,
                     food == "red meat"), mean_percent),
         pull(filter(diet_and_guidelines,
                     food == "red meat"), sex))

    Mood two-sample test of scale

data:  pull(filter(diet_and_guidelines, food == "red meat"), mean_percent) and pull(filter(diet_and_guidelines, food == "red meat"), sex)
Z = NaN, p-value = NA
alternative hypothesis: two.sided

var.test(log10(pull(filter(diet_and_guidelines,
                     food == "red meat"), mean_percent))~
         pull(filter(diet_and_guidelines,
                     food == "red meat"), sex))

    F test to compare two variances

data:  log10(pull(filter(diet_and_guidelines, food == "red meat"), mean_percent)) by pull(filter(diet_and_guidelines, food == "red meat"), sex)
F = 0.96369, num df = 194, denom df = 194, p-value = 0.797
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.7266925 1.2779719
sample estimates:
ratio of variances 
          0.963687 

The p value >.05 for both tests, thus we can conclude that there is not enough evidence to reject the null (no difference in the spread of the distributions), thus we conclude that variance is roughly equal. Great… so we are in pretty good shape for both the \(t\)-test and the linear regression.

So now we can comparing the consumption of red meat by both genders using both a \(t\)-test and a linear regression:

diet_and_guidelines %>%
  filter(food == "red meat") %>%
  lm(log10(mean_percent) ~ sex, data=.) %>%
  plot()

diet_and_guidelines %>%
  filter(food == "red meat") %>%
  lm(log10(mean_percent) ~ sex, data=.) %>%
  summary()

Call:
lm(formula = log10(mean_percent) ~ sex, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0238 -0.2484  0.0052  0.3127  0.6170 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.79887    0.02451  73.382  < 2e-16 ***
sexMale      0.18439    0.03467   5.319 1.77e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3423 on 388 degrees of freedom
Multiple R-squared:  0.06795,   Adjusted R-squared:  0.06555 
F-statistic: 28.29 on 1 and 388 DF,  p-value: 1.77e-07

diet_and_guidelines %>%
  filter(food == "red meat") %>%
  t.test(log10(mean_percent) ~ sex, data= ., var.equal = TRUE)

    Two Sample t-test

data:  log10(mean_percent) by sex
t = -5.3187, df = 388, p-value = 1.77e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2525474 -0.1162266
sample estimates:
mean in group Female   mean in group Male 
            1.798872             1.983259 

Notice how the \(t\) value and the p-value are the same! (Well alomost, the \(t\) value is negative in the t.test() output because the female group is being used as reference group, while the male group is being used as the reference group in lm()). We can fix this using the fct_inorder() function of the forcats package which is all about factors. This function allows us to order the factor by what appears first. In this case “male” appears first, so now our output will match that of the lm() function.

diet_and_guidelines %<>%
  mutate_at(vars(sex), factor)

diet_and_guidelines %>%
  filter(food == "red meat") %>%
  lm(log10(mean_percent) ~ sex, data=.) %>%
  summary()

Call:
lm(formula = log10(mean_percent) ~ sex, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0238 -0.2484  0.0052  0.3127  0.6170 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.79887    0.02451  73.382  < 2e-16 ***
sexMale      0.18439    0.03467   5.319 1.77e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3423 on 388 degrees of freedom
Multiple R-squared:  0.06795,   Adjusted R-squared:  0.06555 
F-statistic: 28.29 on 1 and 388 DF,  p-value: 1.77e-07

diet_and_guidelines %>%
  mutate_at(vars(sex), forcats::fct_inorder)%>%
  filter(food == "red meat") %>%
  t.test(log10(mean_percent) ~ sex, data= ., var.equal = TRUE)

    Two Sample t-test

data:  log10(mean_percent) by sex
t = 5.3187, df = 388, p-value = 1.77e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1162266 0.2525474
sample estimates:
  mean in group Male mean in group Female 
            1.983259             1.798872 

Now they match. Notice that the degrees of freedom also match, both results show 388 degrees of freedom. We are estimating 2 parameters for the linear model the two \(\beta\) coefficients, (the slope and intercept), and for the \(t\)-test we are estimating the means of two groups (males and females). Overall we have two samples (male and female) for each of the 195 countries.

Thus, the overal sample number is: \(n = 195*2 = 390\)

\[df = n - # parameters estimating\] Thus the degrees of freedom can be calculated as: \(df = 390 -2 = 388\)

In the lm() output, the Beta not (\(\beta_{0}\)) which can be interpreted as the intercept or the mean value when sex is not male (so in this case when sex is female) has the same value as the calculated mean of in the t.test() output. The Beta one (\(\beta_{1}\)) can be interpreted as the slope of the regression line or the difference bewteen the means of the two groups. We can also think of this as as the change in \(Y\) with one unit change in \(X\).

If we subtract the means calculated in the t.test() output, we get the value of \(\beta_{1}\) (the slope or the sexMale estimate) of the lm() output!

Mean of males - Mean of females \(1.983259 - 1.798872 =0.184387\)

Cool!

Let’s plot the data to see what is happening:

diet_and_guidelines%>%
  filter(food == "red meat") %>%
  ggplot(aes(x = sex, y = log10(mean_percent))) + 
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(width = .2)+
  stat_summary(fun.y=mean, colour="blue", geom="point", 
              shape=18, size=3)+
  geom_abline(intercept =1.79887 , slope = 0.18439, color = "red") 

# using the slope and intercept from the lm output

We can see that the slope of the linear model is equivalent to the difference bewteen the mean of the two groups, as the linear model line is parallel to the line that could be drawn between the two means.

Also, since we have female and male values each from the same countries, our data is what we called paired. It would be better if we modeled our data that way. First we will use the pivot_wider() function of the dplyr package to get our data in a format that is easier to use for this purpose. To use this function we specify the values that we want to separate into more variables using the values_from argument and we use the names_from argument to specify how we want to separate these other variables. In this case we will make a male and female version of all the other variables specified.

wide_diet <-pivot_wider(diet_and_guidelines, 
                        values_from= c(contains("percent"),
                                       mean, 
                                       upper, 
                                       lower, 
                                       opt_achieved), 
                        names_from= sex)

glimpse(wide_diet)

Observations: 2,925
Variables: 27
$ year_id                   <dbl> 2017, 2017, 2017, 2017, 2017, 2017, 20…
$ location_name             <chr> "China", "China", "China", "China", "C…
$ rei_id                    <dbl> 111, 112, 113, 114, 115, 116, 117, 118…
$ food                      <chr> "fruits", "vegetables", "whole grains"…
$ age_group_name            <chr> "All Available Ages", "All Available A…
$ parameter                 <chr> "continuous", "continuous", "continuou…
$ unit                      <chr> "g/day", "g/day", "g/day", "g/day", "g…
$ direction                 <chr> "low", "low", "low", "low", "low", "hi…
$ lower_optimal             <dbl> 200.0, 290.0, 100.0, 16.0, 350.0, 18.0…
$ optimal                   <dbl> 250.00, 360.00, 125.00, 21.00, 435.00,…
$ upper_optimal             <dbl> 300.0, 430.0, 150.0, 25.0, 520.0, 27.0…
$ unit_optimal              <chr> "g", "g", "g", "g", "g", "g", "g", "g"…
$ food_to_plot              <chr> "fruits", "vegetables", "whole\ngrains…
$ mean_percent_Male         <dbl> 30.1276714, 69.9780163, 20.5969990, 9.…
$ mean_percent_Female       <dbl> 34.3698972, 65.4711643, 19.2030299, 8.…
$ range_percent_Male        <dbl> 37.6595893, 86.8692617, 25.7462488, 12…
$ range_percent_Female      <dbl> 42.9623715, 81.2745488, 24.0037874, 11…
$ percent_over_under_Male   <dbl> 62.34041, 13.13074, 74.25375, 87.46243…
$ percent_over_under_Female <dbl> 57.03763, 18.72545, 75.99621, 88.20590…
$ mean_Male                 <dbl> 7.531918e+01, 2.519209e+02, 2.574625e+…
$ mean_Female               <dbl> 8.592474e+01, 2.356962e+02, 2.400379e+…
$ upper_Male                <dbl> 8.658295e+01, 2.946774e+02, 2.723223e+…
$ upper_Female              <dbl> 9.759191e+01, 2.720864e+02, 2.541951e+…
$ lower_Male                <dbl> 6.543868e+01, 2.153508e+02, 2.426562e+…
$ lower_Female              <dbl> 7.531231e+01, 2.033138e+02, 2.262773e+…
$ opt_achieved_Male         <chr> "No", "No", "No", "No", "No", "No", "Y…
$ opt_achieved_Female       <chr> "No", "No", "No", "No", "No", "No", "Y…

# wide_diet %>%
#   filter(food == "red meat") %>%  
#   lm(log10(mean_percent_Female) - log10(mean_percent_Male) ~1, data= .) %>%
#   summary()

library(lmerTest)

summary(lmer(log10(pull(filter(diet_and_guidelines, food == "red meat"),
                  mean_percent)) ~ 
             pull(filter(diet_and_guidelines, food == "red meat"),
                  sex)+ (1 | pull(filter(diet_and_guidelines, food == "red meat"), location_name))))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
log10(pull(filter(diet_and_guidelines, food == "red meat"), mean_percent)) ~  
    pull(filter(diet_and_guidelines, food == "red meat"), sex) +  
        (1 | pull(filter(diet_and_guidelines, food == "red meat"),  
            location_name))

REML criterion at convergence: -969.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2253 -0.3914 -0.0038  0.3993  3.2527 

Random effects:
 Groups                                                              
 pull(filter(diet_and_guidelines, food == "red meat"), location_name)
 Residual                                                            
 Name        Variance  Std.Dev.
 (Intercept) 1.171e-01 0.342181
             9.363e-05 0.009676
Number of obs: 390, groups:  
pull(filter(diet_and_guidelines, food == "red meat"), location_name), 195

Fixed effects:
                                                                Estimate
(Intercept)                                                    1.799e+00
pull(filter(diet_and_guidelines, food == "red meat"), sex)Male 1.844e-01
                                                               Std. Error
(Intercept)                                                     2.451e-02
pull(filter(diet_and_guidelines, food == "red meat"), sex)Male  9.799e-04
                                                                      df
(Intercept)                                                    1.942e+02
pull(filter(diet_and_guidelines, food == "red meat"), sex)Male 1.940e+02
                                                               t value
(Intercept)                                                      73.38
pull(filter(diet_and_guidelines, food == "red meat"), sex)Male  188.16
                                                               Pr(>|t|)
(Intercept)                                                      <2e-16
pull(filter(diet_and_guidelines, food == "red meat"), sex)Male   <2e-16
                                                                  
(Intercept)                                                    ***
pull(filter(diet_and_guidelines, food == "red meat"), sex)Male ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
p((__,f="ms -0.020

t.test(log10(pull(filter(wide_diet, food == "red meat"),
                  mean_percent_Female)), 
       log10(pull(filter(wide_diet, food == "red meat"),
                  mean_percent_Male)),
       var.equal = TRUE, paired = TRUE)

    Paired t-test

data:  log10(pull(filter(wide_diet, food == "red meat"), mean_percent_Female)) and log10(pull(filter(wide_diet, food == "red meat"), mean_percent_Male))
t = -188.16, df = 194, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1863197 -0.1824542
sample estimates:
mean of the differences 
              -0.184387 

You can see that now our degrees of freedom are 194, which makes sense becuase now we have matched or paired values for the 195 different regions and so \(df = n-1 = 195 -1 = 194\).

Again we get the same \(t\)-value and p-value with both tests.

We can conclude from these tests that we have enough evidence to reject the null hypothesis that there is no difference between the means or that sex has no association or influence on red meat consumption. Therefore, it appears that men consume significantly more red meat than females globally.

Let’s choose another dietary factor to compare between males and females.

diet_and_guidelines_long %>%
  filter(food =="fruits") %>%
  ggplot(aes(x=log10(percent))) +
  geom_histogram() +
  facet_wrap(~ percent_type)

diet_and_guidelines_long %>%
    filter(food =="fruits") %>%
  ggplot(aes(sample = log10(percent))) +
  facet_wrap(~ percent_type)+
  geom_qq() +
  geom_qq_line()

mood.test(pull(filter(diet_and_guidelines,
                     food == "fruits"), mean_percent),
         pull(filter(diet_and_guidelines,
                     food == "fruits"), sex))

    Mood two-sample test of scale

data:  pull(filter(diet_and_guidelines, food == "fruits"), mean_percent) and pull(filter(diet_and_guidelines, food == "fruits"), sex)
Z = 0, p-value = 1
alternative hypothesis: two.sided

var.test(log10(pull(filter(diet_and_guidelines,
                     food == "fruits"), mean_percent))~
         pull(filter(diet_and_guidelines,
                     food == "fruits"), sex))

    F test to compare two variances

data:  log10(pull(filter(diet_and_guidelines, food == "fruits"), mean_percent)) by pull(filter(diet_and_guidelines, food == "fruits"), sex)
F = 1.0025, num df = 194, denom df = 194, p-value = 0.9859
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.7559929 1.3295002
sample estimates:
ratio of variances 
          1.002543 

# wide_diet %>%
#   filter(food == "fruits") %>%  
#   lm(log10(mean_percent_Female) - log10(mean_percent_Male) ~1, data= .) %>%
#   summary()

summary(lmer(log10(pull(filter(diet_and_guidelines, food == "fruits"),
                  mean_percent)) ~ 
             pull(filter(diet_and_guidelines, food == "fruits"),
                  sex)+ (1 | pull(filter(diet_and_guidelines, food == "fruits"), location_name))))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: 
log10(pull(filter(diet_and_guidelines, food == "fruits"), mean_percent)) ~  
    pull(filter(diet_and_guidelines, food == "fruits"), sex) +  
        (1 | pull(filter(diet_and_guidelines, food == "fruits"),  
            location_name))

REML criterion at convergence: -1245.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.4919 -0.3057 -0.0005  0.2940  6.5088 

Random effects:
 Groups                                                            
 pull(filter(diet_and_guidelines, food == "fruits"), location_name)
 Residual                                                          
 Name        Variance  Std.Dev.
 (Intercept) 8.212e-02 0.286573
             3.222e-05 0.005676
Number of obs: 390, groups:  
pull(filter(diet_and_guidelines, food == "fruits"), location_name), 195

Fixed effects:
                                                               Estimate
(Intercept)                                                   1.562e+00
pull(filter(diet_and_guidelines, food == "fruits"), sex)Male -5.570e-02
                                                             Std. Error
(Intercept)                                                   2.053e-02
pull(filter(diet_and_guidelines, food == "fruits"), sex)Male  5.748e-04
                                                                     df
(Intercept)                                                   1.941e+02
pull(filter(diet_and_guidelines, food == "fruits"), sex)Male  1.940e+02
                                                             t value
(Intercept)                                                    76.12
pull(filter(diet_and_guidelines, food == "fruits"), sex)Male  -96.89
                                                             Pr(>|t|)    
(Intercept)                                                    <2e-16 ***
pull(filter(diet_and_guidelines, food == "fruits"), sex)Male   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
p((__,f=="s -0.014

t.test(log10(pull(filter(wide_diet, food == "fruits"),
                  mean_percent_Female)), 
       log10(pull(filter(wide_diet, food == "fruits"),
                  mean_percent_Male)),
       var.equal = TRUE, paired = TRUE)

    Paired t-test

data:  log10(pull(filter(wide_diet, food == "fruits"), mean_percent_Female)) and log10(pull(filter(wide_diet, food == "fruits"), mean_percent_Male))
t = 96.892, df = 194, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.05456356 0.05683104
sample estimates:
mean of the differences 
              0.0556973 

means<-wide_diet %>% 
   filter(food== "fruits") %>% 
   summarise(female =mean(log10(mean_percent_Female)), male = mean(log10(mean_percent_Male)))

means

# A tibble: 1 x 2
  female  male
   <dbl> <dbl>
1   1.56  1.51

Notice here that the mean of the differences value is indeed the same as the coefficient and the \(t\)-value and the \(p\)-value match. The \(\beta_{1}\) coefficient value (slope) matches the difference between the means value calcuated in \(t\)-test. The intercept coefficient \(\beta_{0}\) value of 1.56 (when sex is zero or female) also matches the mean value for the female group when calculated separately.

Influence of age on dietary outcomes

Now we are interested in the influence of age group on dietary consumption.

If we wanted to test the hypothesis that there are anygroup differences, that at least one of the groups is different from the others; we could use an ANOVA test.

It turns out that ANOVAs are also a specialized form of regression.

Let’s take a look at how the consumption of fruits varies by age group.

summary(aov(log10(pull(filter(sep_age_diet_data, food == "fruits"), 
               mean)) ~ 
    pull(filter(sep_age_diet_data, food == "fruits"),
         age_group_name)))

                                                                    Df
pull(filter(sep_age_diet_data, food == "fruits"), age_group_name)   14
Residuals                                                         5865
                                                                  Sum Sq
pull(filter(sep_age_diet_data, food == "fruits"), age_group_name)   10.6
Residuals                                                          481.4
                                                                  Mean Sq
pull(filter(sep_age_diet_data, food == "fruits"), age_group_name)  0.7549
Residuals                                                          0.0821
                                                                  F value
pull(filter(sep_age_diet_data, food == "fruits"), age_group_name)   9.197
Residuals                                                                
                                                                  Pr(>F)
pull(filter(sep_age_diet_data, food == "fruits"), age_group_name) <2e-16
Residuals                                                               
                                                                     
pull(filter(sep_age_diet_data, food == "fruits"), age_group_name) ***
Residuals                                                            
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sep_age_diet_data%>%
  filter(food == "fruits") %>%  
  lm(log10(mean) ~age_group_name, data= .) %>%
  summary()

Call:
lm(formula = log10(mean) ~ age_group_name, data = .)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.02816 -0.15255  0.05033  0.19867  0.65109 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             1.87523    0.01447 129.591  < 2e-16 ***
age_group_name30 to 34  0.01006    0.02046   0.491 0.623188    
age_group_name35 to 39  0.02203    0.02046   1.076 0.281801    
age_group_name40 to 44  0.04245    0.02046   2.074 0.038107 *  
age_group_name45 to 49  0.07538    0.02046   3.683 0.000232 ***
age_group_name50 to 54  0.10493    0.02046   5.128 3.03e-07 ***
age_group_name55 to 59  0.13008    0.02046   6.356 2.22e-10 ***
age_group_name60 to 64  0.13884    0.02046   6.785 1.28e-11 ***
age_group_name65 to 69  0.12605    0.02046   6.159 7.78e-10 ***
age_group_name70 to 74  0.11220    0.02046   5.483 4.36e-08 ***
age_group_name75 to 79  0.09762    0.02046   4.770 1.88e-06 ***
age_group_name80 to 84  0.06638    0.02046   3.244 0.001187 ** 
age_group_name85 to 89  0.06207    0.02046   3.033 0.002432 ** 
age_group_name90 to 94  0.05928    0.02046   2.897 0.003785 ** 
age_group_name95 plus   0.05819    0.02046   2.843 0.004478 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2865 on 5865 degrees of freedom
Multiple R-squared:  0.02148,   Adjusted R-squared:  0.01915 
F-statistic: 9.197 on 14 and 5865 DF,  p-value: < 2.2e-16

Influence of region on dietary outcomes

Data Visualization

Summary

Next, we want to propose our idea for a second case study in the same focus area of:

• Question: How do diets (consumption of various major foods and nutrients) differ by regions (e.g. low-income countries vs not) and by gender around the world? What is the impact of the different diets on mortality or life expectancy?
• Description: We would take some of the main results from this paper (https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(19)30041-8/fulltext){target="_blank"} and highlight them in a case study (similar to before not trying to reproduce the analysis). This paper shows the leading dietary risk factors for mortality included high intake of sodium, low intake of whole grains, and low intake of fruits.
• Why is this important? This analysis demonstrates the need to improve diet across nations and inform implementation of evidence-based dietary interventions.
• Data: This the link Jess Fanzo sent us (https://vizhub.healthdata.org/gbd-compare/){target="_blank"}, but we had to track down a specific set of files that were not included in the link to be able to do this analysis. These files were sent to us in a CSV format and we have been given permission to host them as part of our case study and just cite the original study.
• Major Data Science Objectives:

1. scraping data (gapminder) # probably not anymore
2. loading data from data package
3. wrangling - joining dplyr
4. visualization - ggplot

• Statistics objectives: We are currently torn between demonstrating linear regression or factor analysis. Once we start the case study, it will be easier for us to understand. We’ll keep you posted on that. # to get a ratio of veggie consumption: Veggies$mean/250 *100 #https://www.guru99.com/r-anova-tutorial.html

To do this we need to check if our data violates the assumptions that these statistical test rely on in order to choose the appropriate test to compare groups.